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Summary Report 
of 

THE AIR-BLAST-INDUCED ENVIRONMENT 
WITHIN CIVIL DEFENSE BLAST-SLANTED SHELTERS 

The objective of this work unit was to review and discuss the application 
of hydrodynamic codes in the Office of Civil Defense (OCD) shelter research 
program, and then consider the development of a conceptual model of the flow 
into and through blast-slanted basement shelters. The review selected a 
two-dimensional Eulerian code which numerically simulated inviscid, com¬ 
pressible flow without heat transfer. This form of code was shown to be 
the most suitable form of code for OCD needs on a cost-effectiveness basis. 

The Eulerian equations were discussed, linked to their numerical form, 
and a complete set of numerical equations was appended for the reader's 
convenience. To provide the basis for the evaluation of this approach to 
future OCD problems a brief review of the internal workings of this form 
of code was undertaken. A bibliography was provided which will allow the 
reader to expand the brief review provided by this work unit. 

To test the accuracy of numerical simulation the code was applied to 
several situations which represented conditions that would have to be 
reproduced in the shelter research program. These situations had been simu¬ 
lated experimentally so that a direct comparison between simulation proce¬ 
dures could be undertaken. The first cases reviewed the code's ability to 
reproduce a shock wave profile, and then surveyed the difficulties of 
correctly matching the time frames of rapidly decaying waves between experi¬ 
mental and numerical simulations. An appended section discusses some of 
these difficulties based on a geometric representation of the pressure-time 
profiles. 

The first comparisons of experimental and numerical simulations of 
shelter situations used small scale model shelters. Shelter models used in 
scaled experiments (referenced in text) were reproduced in the numerical 
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simulations to allow direct comparisons. Some difficulties in matching 
pressure—time profiles were noted, but an agreement of better than 90 percent 
between the simulations was noted for the extremely short term transients of 
the filling process. The short term transients were also shown not to affect 
the long term fill history of the shelter. 

Trajectories of points formed by a smoke grid in the same scaled 
experiments were used to test the accuracy of the numerically simulated flow 
field. The difficulty of precisely aligning the space-intensity profiles 
in time between simulations was again found. Procedures were developed to 
improve the comparison techniques. In the evolvement of these procedures it 
was established that the problem of timing comparisons was partly due to 
experimental error. An error analysis of the experiments was undertaken and 
pinpointed some errors in the experimental simulations. In general the 
numerical simulations were concluded to produce flow velocities that were 
about 10 percent higher along the inflow axis and about 5 to 10 percent 
lower in the other regions of the shelter. 

Further evaluations were made against full scale shelters. Data were 
gathered from a full-scale field test using high explosives and from some 
preliminary URS shock tunnel tests. The numerical simulations reproduced 
the average pressure-time profiles within the simulated chambers to within 
90 percent of the experimental value. The evolvement of an error range for 
the experimental data proved to be the limiting factor in developing a more 
precise estimate of variation between simulations. 

The pressure data from URS shock tunnel tests, its evaluation, and data 
and analysis of flow visualization experiments were developed in an appended 
section. The presentation of the flow experiments included a consideration 
for the movement of spheres in the inflow. The flow data generated by the 
experimental and numerical simulations varied on the average less than 10 
percent (within the limits of experimental accuracy). The conclusion reached 
was that the numerical simulation procedures were sufficiently accurate to 
be used in the shelter research program. Some potential problem areas were 
outlined for future study. 
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A section on cost and utilization of numerical simulations was included. 
Comparisons of the costs of numerical versus experimental simulations showed 
numerical simulations to be less costly. A common difficulty with data dis¬ 
play for both simulation techniques was noted. The numerical simulations 
were shown to have the advantage of better output and monitoring control. 

The credibility and reliability of results were considered. The possibility 
of using a simple coarse grid code for general use on a small computer was 
discussed. 

A review of the transient inflow process was undertaken. A review of 
incompressible flow studies provided some insight into the dynamics of jet 
formation. This review was followed by the evolvement of an analytic model 
of the transition flow of a compressible gas into the chamber. The model 
used a quasi-one-dimensional representation of the inflow which was analyti¬ 
cally defined by the method of characteristics. The method showed how 
imbalances in the flow occurred due to the high inflow velocity of the 
adiabatic expansion. The result was in effect a piston which drove the gas 
in the chamber ahead of it. The implications were discussed and verified 
from the experimental data presented in the report. 

A review of the turbulence and jet spreading problem was undertaken and 
a bibliography on it and associated flow problems was developed. Many 
estimation procedures for growth and spread of turbulence were noted. The 
discussion used a simple boundary layer model to represent the spread of the 
flow. The model indicated that the boundary layer and turbulence effects 
were much more pronounced in the scaled experiments than would be expected 
in the prototype situation. The work indicated that the spread and growth 
of turbulence and boundary layers would probably not be of major proportions 
in full scale shelters. 

A general view of associated flow problems was undertaken. The most 
significant of these problems was the effect of the inflow into the shelter 
on the reservoir conditions. The inability of the reservoir to completely 
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replace the effluent mass and energy was shown to reduce the driving pressure 
and also to decrease the driven flow. The comments were related to experi¬ 
mental results. 

A series of numerical simulations were undertaken in support of this 
study. The problem of effective and economic data display arose, hence the 
results of these simulations were verbally described. The simulations 
treated a variety of geometric situations and included multiple room 
geometries. Implications of the observed results were discussed. 

The last section summarized the review of the shelter flow in terms of 
a conceptual model of the flow. The use of a systems representation of the 
flow was used to show that a conceptual model on any general scale was 
probably not a reasonable goal. Some suggestions on how to circumvent this 
difficulty, and what areas of the flow problem might yield the most gains 
in future studies were outlined. 


-4- 




UR51 755-3 


CONTENTS 

Section Page 

ACKNOWLEDGMENTS. ii 

ILLUSTRATIONS. v 

TABLES. vii 

SUMMARY . . ‘. ix 

1 INTRODUCTION. 1-1 

Study Outline. 1-3 

2 COMPUTER FLOW CODE EVALUATION. 2-1 

The Eulerian Method. 2-3 

Difference Equations . 2-9 

3 CODE APPLICATION... 3-1 

Accuracy Considerations . 3-1 

Comparison of Numerical Simulation with Full Scale 

Experimental Results . ... 3-31 

Utilization, Cost, and Display of Numerical Simulations . . 3-36 

Utilization. 3-38 

Cost. 3-41 

Display. 3-43 

4 FLOW PROCESSES ... 4-1 

The Transient Inflow Process . 4-1 

General Comments . 4-33 

5 REVIEW OF NUMERICAL SIMULATIONS.. 5-1 

Chamber Length Variation Simulations . 5-1 

Entrance Location Effects . 5-2 

Multiple Chamber Simulations . 5-4 

Comments. 5-6 

6 DISCUSSION. 6-1 

7 REFERENCES AND BIBLIOGRAPHY . 7-1 

iii 



























URE5 ■ 755-3 


CONTENTS (Continued) 

Appendix Page 

A REPRESENTATIVE FINITE-DIFFERENCE EQUATIONS . A-l 

B A DISCUSSION OF THE GEOMETRIC ASPECTS OF 

PRESSURE MEASUREMENT . B-l 

Introduction . B-l 

Discussion. B-4 

Conclusions. B-18 

C PRESENTATION OF PRELIMINARY URS SHOCK TUNNEL TEST DATA . . C-l 


IV 







LIST OF ILLUSTRATIONS 


UR5 ■ 755-3 

Figure Page 

1 An Example of Two-Dimensional Eulerian Mesh . 2-5 

2 BRL Smoke Grid Model.. 3-2 

3 Pressure Profile Across a Shock Wave...3-3 

4 Velocity Profile Across a Shock Wave . . 3-4 

5 Pressure Profile Across a Reflecting Shock Wave . 3-6 

6 Comparison of Numerical and Model Experimental Simulations , . . 3-9 

7 Particle Velocity Comparison . 3-12 

8 Velocity Time Profiles for Fixed Positions in the BRL Model. . , 3-13 

9 Dimensions of Simulation Field ... . 3-14 

10 Velocity-Time Profiles . ..... . 3-15 

11 Experimental Flow Contours at Time = 215 jUsec.3-22. 

12 Experimental Flow Contours at Time = 405 ^sec.3-23 

13 Experimental Flow Contours at Time = 713 j^sec.3-24 

14 Comparison of Model Configurations.3-29 

15 Comparison of Numerical and Experimental Full Scale 

Pressure Histories . 3-32 

16 A Two Dimensional Simulation of a Free Field Test Case . 3-33 

17 Orientation of Shelter Models with Respect to the Incident 

Shock Waves for Tests in Ref. 43.3-39 

18 Grid for Comparison of Variable-Cell-Size Procedure . 3-46 

19 Grid for Comparison of Variable-Cell-Size Procedure . 3-49 

20 Cumulative Mass Influx as a Function of Time.3-52 


v 




















UR5 ■ 


Figure Page 

21 Cumulative Total Energy Influx as a Function of Time ...... 3-53 

22 Cumulative Potential Energy Influx as a Function of Time .... 3-54 

23 Cumulative Kinetic Energy Influx as a Function of Time . 3-55 

24 Mass Influx as a Function of Computation Intervals ..3-57 

25 Total Energy Influx as a Function of Computation Intervals • • • 3-58 

26 Potential Energy Influx as a Function of Computation Intervals . 3-59 

27 Kinetic Energy Influx as a Function of Computation Intervals . . 3-60 

28 Vena Contracta and Cell Size Geometry ............. 3-62 

29 Curie’s Representations of Early Inflow Process ........ 4-4 

30 Blast Transmission into a Narrow Entrance and the Equivalent 

Process in a Long Entry...4-7 

31 A Simple Diagram of a Quasi-One-Dimensional Process ...... 4-8 

32 Preliminary Wave Diagram of One-Dimensional Model . 4-11 

33 Wave Diagram Approximation for Quasi-One-Dimensional 

Inflow Simulation . 4-14 

34 Schematic of Jet Flow.4-19 

35 Centerline Velocity Distribution . . . 4-21 

36 The Velocity Profiles for a Plane Jet.4-22 

37 Transient Boundary Layer Transformation to Steady Flow 

Analogy..4-27 

38 Approximated Boundary Layer Thickness Parameters for a 

Free-Flow Jet. 4-31 

39 Average Entry Dynamic Pressure ..... . . . 5-7 

40 A Simple Graphic Representation of Flow Field Boundary 

Conditions.6-2 

vi 


















UR5i 


LIST OF TABLES 

Table Pa g e 

1 Flow Vector Comparison at Time = 214 fi sec. 3-18 

2 Flow Vector Comparison at Time = 406 jjisec . 3-19 

3 Flow Vector Comparison at Time = 713 (isec . 3-20 

4 An Outline of Numerical and Experimental 

Simulation Procedures . 3-42 

5 Comparison of Multiply Variable Cells with Non-variable 

Cell Size... 3-47 

6 Comparison of the Effect of Grid Coarseness on the 

Output of Numerical Simulations . 3-50 

7 Comparison of Flow Particle Velocities . 4-10 

8 Flow Field Parameters for Fig. 33. 4-16 


vii 










UIJ 755-3 


NOTATION 

V is the vector velocity 

V is the vector operator 

grad is the vector operator 

div is the vector operator 

p is density 

t is time 

E is total energy 

e is internal energy 

p is pressure 

V is the control volume of an element 

A is surface area of the control volume 

A is the area including the local normal unit vector orienting 

the surface 

a sound speed 

V ratio of specific heats 

S is entropy 
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SUMMARY 

The objective of this work unit was to review and discuss the application 
of hydrodynamic codes in the Office of Civil Defense (OCD) shelter research 
program, and then consider the development of a conceptual model of the flow 
into and through blast-slanted basement shelters. The review selected a 
two-dimensional Eulerian code which numerically simulated inviscid, com¬ 
pressible flow without heat transfer. This form of code was shown to be 
the most suitable form of code for OCD needs on a cost-effectiveness basis. 

The Eulerian equations were discussed, linked to their numerical form, 
and a complete set of numerical equations was appended for the reader f s 
convenience. To provide the basis for the evaluation of this approach to 
future OCD problems a brief review of the internal workings of this form 
of code was undertaken. A bibliography was provided which will allow the 
reader to expand the brief review provided by this work unit. 

To test the accuracy of numerical simulation the code was applied to 
several situations which represented conditions that would have to be 
reproduced in the shelter research program. These situations had been simu¬ 
lated experimentally so that a direct comparison between simulation procedures 
could be undertaken. The first cases reviewed the code T s ability to reproduce 
a shock wave profile, and then surveyed the difficulties of correctly matching 
the time frames of rapidly decaying waves between experimental and numerical 
simulations. An appended section discusses some of these difficulties based 
on a geometric representation of the pressure—time profiles. 

The first comparisons of experimental and numerical simulations of 
shelter situations used small scale model shelters. Shelter models used in 
scaled experiments (referenced in text) were reproduced in the numerical simu¬ 
lations to allow direct comparisons. Some difficulties in matching pressure- 
time profiles were noted, but an agreement of better than 90 percent between 
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the simulations was noted for the extremely short term transients of the 
filling process. The short term transients were also shown not to affect 
the long term fill history of the shelter. 

Trajectories of points formed by a smoke grid in the same scaled 
experiments were used to test the accuracy of the numerically simulated flow 
field. The difficulty of precisely aligning the space-intensity profiles 
in time between simulations was again found. Procedures were developed to 
improve the comparison techniques. In the evolvement of these procedures it 
was established that the problem of timing comparisons was partly due to 
experimental error. An error analysis of the experiments was undertaken and 
pinpointed some errors in the experimental simulations. In general the numer¬ 
ical simulations were concluded to produce flow velocities that were about 
10 percent higher along the inflow axis and about 5 to 10 percent lower in 
the other regions of the shelter. 

Further evaluations were made against full scale shelters. Data were 
gathered from a full-scale field test using high explosives and from some 
preliminary URS shock tunnel tests. The numerical simulations reproduced the 
average pressure-time profiles within the simulated chambers to within 90 
percent of the experimental value. The evolvement of an error range for the 
experimental data proved to be the limiting factor in developing a more 
precise estimate of variation between simulations. 

The pressure data from URS shock tunnel tests, its evaluation, and data 
and analysis of flow visualization experiments were developed in an appended 
section. The presentation of the flow experiments included a consideration 
for the movement of spheres in the inflow. The flow data generated by the 
experimental and numerical simulations varied on the average less than 10 
percent (within the limits of experimental accuracy). The conclusion reached 
was that the numerical simulation procedures were sufficiently accurate to 
be used in the shelter research program. Some potential problem areas were 
outlined for future study. 
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A section on cost and utilization of numerical simulations was included. 
Comparisons of the costs of numerical versus experimental simulations showed 
numerical simulations to be less costly. A common difficulty with data dis¬ 
play for both simulation techniques was noted. The numerical simulations 
were shown to have the advantage of better output and monitoring control. 

The credibility and reliability of results were considered. The possibility 
of using a simple coarse grid code for general use on a small computer was 
discussed. 

A review of the transient inflow process was undertaken. A review of 
incompressible flow studies provided some insight into the dynamics of jet 
formation. This review was followed by the evolvement of an analytic model 
of the transition flow of a compressible gas into the chamber. The model 
used a quasi-one-dimensional representation of the inflow which was analyti¬ 
cally defined by the method of characteristics. The method showed how 
imbalances in the flow occurred due to the high inflow velocity of the 
adiabatic expansion. The result was in effect a piston which drove the gas 
in the chamber ahead of it. The implications were discussed and verified 
from the experimental data presented in the report. 

A review of the turbulence and jet spreading problem was undertaken and 
a bibliography on it and associated flow problems was developed. Many 
estimation procedures for growth and spread of turbulence were noted. The 
discussion used a simple boundary layer model to represent the spread of the 
flow. The model indicated that the boundary layer and turbulence effects 
were much more pronounced in the scaled experiments than would be expected 
in the prototype situation. The work indicated that the spread and growth 
of turbulence and boundary layers would probably not be of major proportions 
in full scale shelters. 

A general view of associated flow problems was undertaken. The most 
significant of these problems was the effect of the inflow into the shelter 
on the reservoir conditions. The inability of the reservoir to completely 
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replace the effluent mass and energy was shown to reduce the driving pressure 
and also to decrease the driven flow. The comments were related to experi¬ 
mental results. 

A series of numerical simulations were undertaken in support of this 
study. The problem of effective and economic data display arose, hence the 
results of these simulations were verbally described. The simulations treated 
a variety of geometric situations and included multiple room geometries. 
Implications of the observed results were discussed. 

The last section summarized the review of the shelter flow in terms of 
a conceptual model of the flow. The use of a systems representation of the 
flow was used to show that a conceptual model on any general scale was 
probably not a reasonable goal. Some suggestions on how to circumvent this 
difficulty, and what areas of the flow problem might yield the most gains 
in future studies were outlined. 
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Section 1 
INTRODUCTION 


One option being investigated for the future design and utilization of 
civil defense shelters is to keep the shelter open during an attack situation. 
The open entrances allow the blast environment to be transmitted into the shel¬ 
ter in the form of potentially hazardous high-velocity flows. These flows are 
of interest to the designer because of the structural loads that they form and 
their potential for causing injury to the shelter occupants. 

The damage potential of the inflow into shelters for incident over¬ 
pressures below 15 psi was documented in Refs. 1 and 2. These studies further 
showed that the induced flow field could be the major cause of human injury. 
The basis for the work in Ref. 2 was two earlier studies (Refs. 3 and 4) which 
first continued previous work on the prediction methods for chamber filling, 
and then the flow field within the chambers. These definitions illustrated 
the intensity of the inflow, provided some upper bounds on the magnitude of 
the flow-field parameters, and established where additional work was required 
to better define the inflows. 

The second generation of experimental and analytical studies (Refs. 2 
and 5) better defined the flow field and evaluated the effectiveness of flow- 
limiting devices. Within these OCD work units, some better estimates of 
inflow parameters were determined. Some of the previous problem areas did 
recur, namely, the description of the flow field during the early part of 
the filling process, the establishment of a model for the spreading of the 
inflow, and the development of circulating flows within the shelters. 

To provide the shelter designer with more accurate information, all of 
these studies had to be linked together. The concept of the open-door, blast- 
slanted basement shelter (Ref. 6) requires the flow to be defined in greater 
detail since the feasibility of the concept is dependent upon the casualty 
expectations from all causes including the high velocity flows. The preceding 
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studies bounded the problems but did not provide the complete definition of 
the flows required to ascertain the feasibility of the concept. This study 
undertook to meet this specific need of the Office of Civil Defense planners 
by providing improved understanding of the flow as well as linking all the 
existing information into a conceptual model of engineering accuracy. 

The study of these flows is complicated by the wide range of geometries 
encountered in basements, and by the transiencies of the flows. The component 
flows which make up the total flow field also change and interact. This com¬ 
plexity precludes an explicit analytic solution. The problem has been 
classically approached either through experimental studies, through analytical 
descriptions of the component flows (usually using quasi-steady-state 
approximations), or through numerical simulations. The final engineering 
solution to these qualitative and quantitative descriptions of the flows will 
be made through a synthesis of all three approaches. 

Experimental studies will provide the data against which the analyti¬ 
cally based techniques must prove themselves. The limitations of experimental 
procedures in the greatly space-time-dependent flows preclude an empirical 
solution, if any reasonable cost-effectiveness-benefit ratios are to be main¬ 
tained. The variation in potential geometries would increase the necessary 
testing to unreasonable levels. 

The use of numerical simulations poses similar problems, i.e., an 
excessive number of simulations is required to provide an empirical solution 
from the data. The numerical simulations do provide considerably more infor¬ 
mation on the space-time variable flows, in terms of both the flow and 
thermodynamic variables, at a lower cost per unit of information than the 
experimental methods. Even so, the numerical simulations will still not 
provide the overall solution required. 

The use of the individual flow types presents analytic expressions for 
the flows. These can be approached as quasi-steady-state solutions and then 
pieced together. The quasi-steady solutions are not quite precise in the 
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formative stages of the flows; hence analytic models to describe the vari 
ations must be developed. 


The numerical simulations are useful in this context because they pro¬ 
vide a review of the formative flows and describe how they relate to form 
the entire flow field. The experimental data can then be used to review and 
check the accuracy of the simulations and of the conceptual flow model. The 
last stage is to relate quantitative prediction schemes to the qualitative 
model of the flow. The work to be described in this report used the above 
approach as described by the following statement of objectives and scope from 
the OCD work unit outline: 

OBJECTIVE 

The objective of the proposed work is to define (relative to 
the potential hazard to shelter occupants) the characteristics 
and intensity of the flows formed within civil defense base¬ 
ment shelters by the transmission of an air blast environment 
into the shelter. 

SCOPE 

The study will develop a conceptual model for the description 
of air blast phenomena inside an open shelter. The form of 
the conceptual model will allow its use in defining engineer¬ 
ing equations and empirical relationships to describe the 
environment within the shelter during blast loading. In 
addition to analytical inputs, the study will use hydro- 
dynamic computer code calculations to develop the con¬ 
ceptual model. As part of the work unit task, a brief review 
of existing codes will be made and a suitable code will be 
selected to perform the computer calculations and evaluated. 


STUDY OUTLINE 

The initial emphasis of the work unit was to consider the selection of 
the computer code and provide a brief description of code formulation and 
function. The selected code work was then applied to problems delineated in 
past work to determine the accuracy of numerical simulations relative to the 
shelter flows being studied. Concurrent with this effort, analytical studies 
were undertaken to begin developing the conceptual model. The computer code 
was used to support these studies as needed. 
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Because of funding constraints, the study had to rely strongly on 
development of bounds on the magnitudes of the flows. Primary emphasis was 
placed on attempts to develop a conceptual model and the computer code as out¬ 
lined in the scope. The studies of particular facets of the flow were 
limited in scope and depth to a large degree by the funds available after the 
primary objectives had been met. These limitations affected the breadth and 
depth of the conceptual model, though the model presented is valid for 
engineering estimates. 

The attempt to develop a conceptual model was undertaken as much as 
possible from a total systems viewpoint and, therefore, attempted to consider 
parameters such as blast wave-structure orientation, architectural/engineering 
constraints, and shelter management procedures. The inclusion of these other 
parameters forced the flow model into the proper perspective, that of being 
but one part of a larger system which controlled the flow f s boundary 
conditions. The application of this technique was somewhat restricted under 
the above-mentioned funding constraint, but its use was felt to be necessary 
in the development of the most practical and realistic conceptual model of 
the flow in the shelter. 
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Section 2 

COMPUTER FLOW CODE EVALUATION 


In order to select a code to simulate the transmission of the air blast 
environment into the basement shelters, the constraints developed by flow- 
shelter environment must be considered. The incident and transmission flows 
are compressible flows with shock waves. The range of incident shock over¬ 
pressure for open blast-slanted shelter design extends to about 15 psi ; hence, 
heat transfer can be neglected. The flows into the shelters form jets, and 
turbulent momentum transport does occur in the transmission flow. 

The shelter doors extend nearly the full height of the room, making flow 
in the shelter approximately two-dimensional. The various entrance confi¬ 
gurations are in reality three-dimensional; but, by careful selection of 
geometries, the entries can be represented in two dimensions. The viscous 
interaction of the flow with the bounding walls is neglected. The lack of 
functional three-dimensional computer codes precludes their application to 
shelter research; therefore, the disadvantages of the two-dimensional approach 
will have to be accepted and evaluated to determine the range of deviations 
they create. 

The techniques available to compute the properties of compressible flow 
into chamber, including shock waves, for two-dimensional space-time 
coordinates are limited. One of the most direct methods uses numerical 
step-wise procedures to compute the flows through zones or cells in the flow 
field by means of finite-difference representations of the equations of motion 
in either a Lagrangian or Eulerian mesh. The zones of the Lagrangian system 
move and distort with the flow. This feature allows particles to be followed 
in space and time, but significant inaccuracies are developed when large 
fluid distortions occur. The Eulerian mesh is a fixed grid in space and time 
through which this fluid flows. This representation allows greater fluid 
distortion without loss of accuracy, but it cannot distinguish between gases 
and particles within the flow field. Some codes, such as the PIC code, 
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combine the two systems, using the best features of each; but they require 
a large computer memory and have long running times. 

The strongly space-time-dependent flows found in the basement shelters 
introduce considerable fluid distortion. The Lagrangian system of 
representation is therefore not a good coordinate system for these problems. 
The disadvantage entailed by the Eulerian system’s not being able to 
distinguish between gases within the flow is not significant. Air is the 
sole gas involved, and since the shock strengths are relatively low, entropy 
increases are minor and differences across contact surfaces are minimized. 
Knowledge of particle location is unnecessary in shelter research; hence, use 
can be made of a strict Eulerian system, with its reduced computer memory and 
running time requirements. 

URS had previously undertaken a similar review process of available 
codes and had established that the Eulerian system would best meet the 
requirements of general two-dimensional shock-wave interaction problems. 

The selection process not only considered the preceding factors, but also 
took into account the machine time and storage required to execute any given 
problem. In that these parameters represent costs to be incurred by the 
project, it was necessary to consider cost-effectiveness parameters in the 
selection process. 

The specific code URS selected was the Fluid-In Cell (FLIC) technique 
developed by Gentry, Martin, and Daly (Ref. 7) of Los Alamos Scientific 
Laboratory to compute two-dimensional time-dependent compressible flows. The 
FLIC code uses a Eulerian mesh of cells in which improved (with respect to 
earlier Eulerian codes) energy, mass, and momentum transport calculations are 
used to reduce the computer memory requirements and to increase the 
computational speed. The flexibility of the FLIC code is enhanced by 
allowing boundary and initial conditions to be input for each specific 
problem without code modifications. 

The development of the FORTRAN IV version of the FLIC code was undertaken 
by URS as an in-house project during recent years. It had been utilized in a 
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previous OCD work unit and shown to provide considerable information within 
engineering accuracy (approximately accurate to i0.9 of the actual value). 
Within the spectrum of Eulerian codes there seems little significant varia¬ 
tion as to their properties. The review of the available Eulerian codes 
undertaken by this study reconsidered the selection process relative to OCD 
shelter problem needs and also included a review of more recently developed 
hydrodynamic computer codes. 

The review did not find another method that would provide any 
significant advantages over the FLIC technique. The evaluation also weighed 
the cost of developing a new code against the cost of the already operational 
FLIC code for use in this project. Based on the review and the cost 
constraints, the FLIC code was selected for use in order to conserve funds. 

A description of general code capabilities and availability relative to 
Office of Civil Defense needs will be presented to shorten the review process 
for future problems. 

THE EULERIAN METHOD 

The problem of simulating a general compressible flow is complicated by 
several factors, as pointed out by J. Von Neuman and R. D. Richtmyer in their 
pioneer work on hydrodynamic codes (Refs. 8 and 9): 

• The presence of shock waves across which the thermodynamic 
variables are discontinuous 

• The difficulty of relating boundary conditions across 
these discontinuities 

• The space-time variation of the shock waves and of the 
flow field itself 

• The motion of the waves and the flow-field composition 
are not known in advance; hence the boundary conditions 
on the shock waves are space-time dependent 

Von Neuman and Richtmyer proposed a method of introducing an artificial 
viscosity term into the difference equations of motion to smooth out the 


2-3 




UR5 ■ 755-3 


discontinuity and make the equations stable. This procedure is the basis 
upon which succeeding codes were developed. 

By means of the smoothing factor, the equations of motion can be written 
in a general form. A given flow field can be divided into sections called 
zones or cells and a set of initial and boundary conditions specified. The 
generalized equations can then be applied sequentially and the flow field 
computed in monotonically increasing time. The modern high-speed digital 
computer reduces to practical limits the computation times and intricate 
bookkeeping required. 

The following discussion of the Eulerian computation sequence for a 
planar two-dimensional isentropic compressible flow will use a mesh of equal 
spacing. The programmer defines the mesh size, the geometry of the field, 
and the boundary and initial conditions. The boundary conditions include the 
boundaries of the field as well as the boundaries of objects within the field 
(see Fig. 1). 

A mesh generation and bookkeeping process must be evolved to control 
the computational procedure. Once the mesh is established and the field 
boundary defined, the process of defining objects within the field must be 
undertaken. The bookkeeping must insure flow is stopped at all solid 
boundaries and must keep the order of the cells. When this is done, the 
computation of the flows can be undertaken. 

Euler’s equations of motion, which form the basis for the computational 
sequence to be described, must fulfill the following set of conditions 
everywhere in the flow field (Ref. 10) . Every particle in the fluid must obey 
Newton’s Second Law at all times. The boundary conditions must be fulfilled 
at the boundaries of each cell. The continuity equations must be satisfied 
at all points in the fluid except at singular points. The importance of 
these conditions will be noted in the later description of the computational 
phases of the Eulerian codes. 
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Fig. 1. An Example of Two-Dimensional Eulerian Mesh 
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The basic equations of motion are those stemming from the equation of 
state and the conservation of momentum, energy, and mass. There are five 
fundamental parameters described as a function of two spatial coordinates at 
any instant of time (Ref. 11). 


5p 

3t 

5v 

at 

as 

at 


- + div OV = 0 

+ V • Vm + — grad p = 0 
0 

+ v • grad s = 0 


P = p(p,s) 


Mass (1) 

Momentum (2) 

Energy (3) 
State (4) 


where S is the entropy, and E, the total energy per unit mass, would be defined 

2 

by E = e + 1/2 (v) . The use of the d f Alembert-Euler acceleration formula 
(Ref. 12) in Eq. (2) should be noted, i.e., 

acceleration = ~ + V ’ Vv (5) 


where ~ is the local acceleration and v * Vv is the convective acceleration, 
ot 

This representation shows both the spatial and temporal variations of 
acceleration with the flow field and is useful for interpreting the flow 
field characteristics to be discussed later in the report. 


There is a more convenient method than the spatial presentation above 
for developing the difference equations for the Eulerian codes. The 
corresponding hydrodynamic equations in Eulerian coordinates for compressible 
isentropic flow in the integral form as defined by Rich (Ref. 13) are: 


/ ft p dV = -£ pv • 

V 

It / pV dV = -£p dA 

V 


(conservation of mass) (6) 


(conservation of momentum) (7) 
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f e dV = - pv • dA (conservation of energy) (8) 

dt J Jk 

v 

p = p(p,E) (equation of state) (9) 

The integrals are evaluated over the volume V of an element of surface A. 
These equations are commonly referred to as the conservation or divergence 
equations. 

A brief review of the properties of the equations will be presented 
since they form the basis for the computational scheme used. Equation (6), 
the conservation of mass (or in this form sometimes called the "Reynolds 
transport theorem,” states that the time rate at which mass increases within 
the volume is equal to the net rate of mass flux across the surface of the 
volume. In terms of the Eulerian calculation grid, each cell becomes a 
volume of unit depth with the surface which must meet the above conservation 
of mass condition. 

Equation (7), the integral expression for the conservation of momentum, 
indicates the condition that at any given time the rate of change of momentum 
in the control volume must equal the resulting force on the volume. The 
resulting force is the net stress or pressure acting on the surface area. 

The last conservation condition, the energy equation as given by Eq. (8), 
expresses the similar conditions for energy that the Reynolds transport 
theorem expressed for mass. The conservation of energy for a control volume 
at any given instant requires that rate of accumulation of energy in the 
form of heat is equal to the rate of energy increase within the control 
volume and the rate work is done by the boundaries on the surroundings. The 
direct relationship between the integral form and its control volume 
representation to the mesh of cells is helpful in understanding the 
computational basis of the Eulerian codes. 

The preceding basic equations are transformed into difference equations 
which correspond to the flow field grid. The transformation occurs in 
relation to a prescribed calculation sequence, which will be described in the 
following paragraphs. The relationship of the sequence to the conditions for 
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an Eulerian system of equations described above should be noted. The 
sequence will be considered, but the difference equations will not be 
presented in the body of this report. For the reader's convenience a full 
set of Eulerian Difference Equations for two-dimensional flow, from the 
work of Rich (Ref. 13), is presented in Appendix A. A brief discussion of 
the differencing and stability criteria for the equations will be presented 
after the computational sequence is considered. 

The first phase of the calculation procedure fixes the fluid within the 
cell, i.e., there is no movement of fluid across the boundaries of the cell. 

The conditions of the preceding step are applied to the fluid element as 
boundary and initial conditions. By applying the equations of motion and 
state, a new set of values (termed intermediate values and defined by the tilde 
symbol) for the velocities and internal energy within the cell can be calculated. 
The procedure is applied to every cell in the mesh. 

The second phase of the calculation procedure allows flow to occur 
across cell boundaries based on the intermediate values calculated in phase 
one. The mass and energy flow out of each cell is tabulated. During this 
phase energy, mass, and momentum are conserved, but no change in the flow 
or state variables is calculated. 

The third phase of the computation sequence sums the mass and energy 
interchange carried out in phase two. The first step of the calculation 
determines a new value of density for each cell in the mesh. By applying 
conservation of energy and momentum to the new flow and state variables, 
values of the flow velocity and specific energy in the cell can be calculated. 

The sequence of operations must occur at all cells within the mesh that 
have open boundaries on all sides. For cells which do not have four open 
sides, a series of new conditions exists. One of the two types of these 
cells uses completely closed boundaries, i.e., the cell is an obstacle 
through which no flow occurs. The other type of cell is a partial cell, in 
which a portion of it allows flow to occur and the remainder is a closed 
cell. These three types of cells establish the flow field and its boundary 
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conditions. The control portion of the code program insures that these 
boundary conditions are adhered to during the computational sequences. The 
only remaining information required to begin the computational sequence is 
the initial conditions of the flow field. Use of the code thus permits one 
to ascertain (in monotonically increasing temporal sequence) the flow field 
properties. 

To better understand the validity and meaning of the computed results, 
the following items will be considered. 

• Differencing techniques 

• Expansions and errors 

• Stability 

• Internal dissipation 

• Cumulative errors 

The first three items are problems which must be controlled if the code is 
to work at all. The latter two items may exist at undesired levels even if 
a stable difference scheme exists. 

DIFFERENCE EQUATIONS 

There are many different approaches to the difference equations which 
can be applied to the general hydrodynamic internal equations. The various 
differencing techniques and their advantages and shortcomings can be found 
in discussions in Refs. 13 through 35 and will not be considered here. The 
following discussion will, however, review those techniques commonly used 
in the two-dimensional hydrodynamic codes and, in particular, those of the 
FLIC code employed in this study. The review will be based on the work in 
Refs. 13 through 17 and will present only that material needed to provide a 
background for the later discussion of code applicability to the OCD type of 
aerodynamic problem. 

It was shown earlier that the conservation equations could easily be 
applied to the individual cells of a finite-difference mesh. The equations 
presented are used in the hydrodynamic codes in their equivalent numerical 
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difference expressions which are obtained by numerical forms at a given 
time. The numerical forms of the equations are then applied with respect 
to the calculational sequence discussed earlier. 

The final phase of the calculation sequence closes the cell boundaries 
to inflow; hence the mass and average density of the fluid in the cell must 
remain constant. The density term can be removed from the time derivative 
in the control volume integration of the conservation of momentum and energy 
equations. By using simple linear averages for the pressure and velocities 
at the boundaries between cells, difference expressions for the time 
derivatives can be written and the resultant expressions for the intermediate 
velocities and specific energies can be written. 


The FLIC code used in this review was developed to allow the use of 
either planar or cylindrical coordinates. To allow for the divergence of 
the cylindrical field in the radial direction, the FLIC equations resemble 
those of the Eulerian portions of the PIC codes (Ref. 14). The differencing 
technique selected for FLIC also varies slightly from the equations presented 
in Appendix A, particularly for the energy equation which follows the 
Zip-type differencing format described by Harlow (Ref. 14). 


The authors of the FLIC code point out that the difference between the 

FLIC format and that presented in Appendix A is the use of forward—backward 

time averages for velocity values, i.e., u^ 1 . = i (u 11 + it 11 ) . The use of 

2 xj lj 

this differencing format limits the occurrence of negative internal energies 
in the computational procedure and also reduces the tendency to destroy 
entropy. The occurrence of negative internal energies results from 
computational instabilities and fluctuations, which tend to decrease the 
internal energies in individual cells. Harlow (ftef. 14) presents the condition 



£ 2(y-l) u 


<h * h) 


( 10 ) 


as the criterion for limiting the occurrence of negative internal energies 
for a polytropic equation of state, p = (y-1) pe. 
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In comparing the equations presented in Appendix A with those of 
Gentry, Martin, and Daly in their presentation of the FLIC code (Ref. 7), the 
appearance of the artificial viscosity term in the pressure term of the 
equations can be noted. The use of the artificial viscosity term enhances 
the stability of the equations in the presence of shock waves in the flow 
field. The value of the term is controlled by the conditions at the 
boundary interfaces and by an arbitrary constant. A more thorough discussion 
of the application procedure can be found in Ref. 7. 

The form of the artificial viscosity, or damping factor, in the 
equation varies quite widely. The methods that can be used are discussed 
and compared by Emery in considerable detail in Ref. 18. The von Neuman- 
Richtmyer form of the damping factor need not be used. For example, the 
equations presented in Appendix A do not include the use of the specific 
external damping term, but rather use a term inherent to the Taylor series 
expansion terms. The FLIC code on the other hand does use the von Neuman- 
Richtmyer artificial viscosity term as modified by Landshoff (Ref. 15) . 

The artificial viscosity term normally is needed in regions where the 
gradients of flow field properties are steep. The strength of the 
dissipative mechanism is a function of the velocity and the mesh size. 
However, neither of these conditions is an absolute guide. For example, 
the effect of artificial damping must be reduced in rarefaction regions and 
enhanced in stagnation regions. To meet all these conditions, various 
dissipative terms are used. The application of the proper term is 
accomplished by the computational controls through a series of test 
conditions. 

The use of the artificial viscosity term enhances the stability of the 
difference equations, but it does not fully eliminate the oscillations behind 
the shock. Landshoff 1 s method does provide fairly accurate representation of 
steep gradients. Relative to the needs of the OCD studies of chamber 
filling, the steepness of the gradients is a useful feature, and the small 
oscillations are a minor disturbance or perturbation in the presentation of 
data. 


2-11 



UR5i 


755-3 


The use of this procedure does introduce the bound on the computational 
time step discussed above, i.e., the maximum allowable time step is a 
function of the cell size and the maximum flow velocity. To insure 
stability, the following further conditions are necessary (Ref. 7): 


Ax 


At < 0.4 


(ID 


cAt 

Ax 


X . 
min 


1 

where \ m = Min[l < B, - YB + (y 2 B 2 + 4y) 2 , YB/(Y-1)] 


( 12 ) 

(13) 


and B is a constant used to determine the magnitude of the viscous pressure 
term and c is sound speed. The inequality expressed in Eq. (11) was found to 
cause some minor difficulties in some instances, and a 0.3 value on the right- 
hand side was found to be a more reliable condition. This equation sets the 
condition that a mass point in a cell cannot flow more than 0.3 AX in one time 
step (At), thus insuring that a cell cannot empty itself. These conditions are 
of practical interest for shelter research since they provide one of the main 
limits on the computational speed of the codes. 


Variation of the FLIC equations, as presented in Ref. 7, from the 
equations presented in Appendix A occurs in phase 2 of the computational 
sequence. The main difference is found in the determination of the flow- 
parameter values at the boundaries of the cells. The FLIC code uses the 
donor-cell method of differencing, in which the transport properties are 
proportional to the properties of the donor cell. The authors indicate an 
increased stability when used in conjunction with the artificial viscosity 
terms. The preceding factors are also used in the last computational phase. 

The preceding discussion has shown dissimilarities exist between the 
various types of differencing techniques. Of primary interest to the OCD 
shelter program is the effect these dissimilarities have on the computational 
accuracy of the blast-induced flows into and through the shelters. Further 
review of the numerical equations will not be undertaken, only a brief 
outline of their related characteristics. Similarly, the development of 
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equations for boundary cells and partial cells will not be reviewed. The 
attached bibliography provides a full description for the reader who wants 
a more thorough coverage and more recent innovations. 

The Eulerian difference equations, as used in the more common numerical 
techniques, have been shown to conserve mass, energy, and momentum (Ref. 13). 
Non-conservative conditions could result due to computational errors, which 
may simply be cumulative errors in the large number of calculations made. 
These errors are easily checked by keeping a running tabulation of the total 
mass, energy, and momentum in the system. If the tabulation remains 
consistent with the accuracy of the computing machine, one check on the 
computational accuracy has been made. 

Even if exact conservation is maintained, errors can exist in the 
numerical solutions. Local perturbations, mathematical instabilities, 
spreading of discontinuities, expansion and internal dissipation, 
mathematical errors, and exclusion of some non-ideal effects can cause 
temporal deviations from the simulation's real world counterparts. It is 
this type of variation that is of most significance in the consideration of 
shelter flows. 

In the initial discussion of the FLIC equations (Ref. 7), the authors 
discussed the errors in the finite-difference equations by expanding them 
in a Taylor series. In reviewing the expansions, the error terms were 
considered to cause artificial diffusion of mass, momentum, and energy where 
sharp gradients in flow properties occur. The physical effect is to cause 
the gradients to spread out spatially. 

The distention of the gradients is a function of the cell size. For 
example, the thickness of a shock wave is usually spread over 3 to 5 cells 
by the dissipative factors. The spread factor for the von Neuman Richtmyer 
artificial viscosity is given by [8B/(Y + l)] 2 , where B is a critical 
constant in the damping, ranging in value between 0.2 and 0.4. To get the 
best resolution, the cell mesh must be kept small. The smaller the spatial 
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or time differences, the smaller the error in the finite-difference 
approximations. By following this argument, infinitely small cells could be 
used, but a limit in available computer memory and computational time is 
soon reached. 

Another source of error is violation of the principle of Galilean 
invariance (the relative motions of particles are the same irrespective of 
the frame of reference). This principle is particularly important when a 
shock wave moves into a flow field, as opposed to a stationary gas. At 
equal Mach numbers, the profiles of the flow parameters across the shock 
wave should be identical. Reference 7 has shown Galilean invariance is 
violated for a shock wave propagating into a fluid moving at one-half the 
shock wave’s velocity. A greater spatial diffusion occurred due to increased 
velocity relative to the computing mesh. In the general case, the diffusion 
is a function of the value of the parameters involved (Mach number, cell 
size, gas properties, etc.). In shelter studies, these variations must be 
considered because strong gradients do exist within counterflows. 

Other regions in which the diffusion terms can cause deviations are 
stagnation and rarefaction regions. In stagnation regions, there is no flow; 
hence the artificial viscosity which is proportional to flow velocities is of 
negligible strength. It has been shown (Ref. 7) that if the conditions expressed 
by Eqs. (12) and (13) are met, then a stable solution is possible. The 
condition is generally the same for all small-magnitude flows. Similarly, 
adjustments for the artificial viscosity must be used in rarefaction regions. 
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Section 3 
CODE APPLICATION 

The discussion in Section 2 reviewed some of the general background, 
limitations, and structure of the Eulerian codes. The usefulness of the 
codes in shelter applications depends on a number of factors: the ability 
to accurately reproduce the blast-induced flow environment in the shelter, 
the cost of reproducing the environment, and the ability to display and uti¬ 
lize the generated data. The accuracy of the code is the most important of 
the factors, because if it does not exist, cost and utilization are only 
academically interesting. Hence the discussion of application will begin 
with a consideration of accuracy. 

ACCURACY CONSIDERATIONS 

Since accuracy is related to a comparison of the numerical simulations 
with prototype conditions, a group of experiments was numerically simulated. 
The first series of these simulation-to-experiment comparisons was made using 
scaled data gathered at the Ballistic Research Laboratories (Ref. 4) for OCD, 
These experimental data relate to the model shown in Fig. 2. Both pressure- 
time histories and smoke grids to measure particle velocities were recorded 
for the model. 

The evaluation will begin by reviewing the code's simulation of a 
classical step shock wave. The pressure and velocity spatial profiles for 
a 1/8-in. mesh size are shown in Figs. 3 and 4. Since the profile is a 
function of the number of cells, finer mesh would provide a correspondingly 
reduced spatial profile. The overshoot and damping characteristics of the 
code can also be observed. 

The steepness of the change of pressure and velocity variables in these 
profiles approaches 1/4 in. for the cells shown. This rise is directly 
comparable to the rise-time associated with the pressure sensors used by 
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Fig. 3. Pressure Profile Across a Shock Wave 
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BRL in the experiments. The overshoot added another 1/8 in. to the profiles 
before the "step" condition was reached. The earlier work (Ref. 2) in this 
area showed this to be of negligible importance to the shelter studies 
because short-term pressure transients did not play a major role in damage 
predictions. 

Another advantage gained by the steep profiles is that they allow the 
use of coarser grids in the simulation procedure without degrading resolution. 
During the period between these reports, the FLIC code was modified to 
accommodate variable cell sizes to increase local resolution without increase 
ing memory requirements (to be discussed more fully in a later section). 

The ability to use finer grids where resolution was required put still less 
of a premium on profile steepness. 

A stability analysis was undertaken by independently varying the coeffi¬ 
cients of the artificial viscosity terms, the cell size, and the computational 
time step increments. From this analysis a smoother, though slightly less 
steep (about 3 to 4 cells), rise was selected. A pressure—distance—time 
profile series is presented in Fig. 5 for a regular reflection process in a 
mesh incorporating two cell sizes. 

The test shown in Fig. 5 is also useful in evaluating the effect of 
mesh size changes, the Galilean invariance, the code properties at reflective 
boundaries, and the stability of the code in low-flow regions. With the 
exception of the Galilean invariance, the test was favorable. The spatial 
spreading of the shock front, although present, is not of unusual proportions 
in comparison to similar results presented in the literature. 

Although the preceding provides some insight into the general accuracy 
of the codes, there are specific questions that must be considered concerning 
the accuracy of the codes relative to the reproduction of flows in basement 
shelters. Attention is being directed to the main damage mechanism (the flow) 
rather than to the pressure profiles of transient shock waves. One source of 
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information available from the BRL experiments is the individual pressure 
histories within the chamber, which provide an indirect measure of one of 
the flow-field parameters. 

These pressure histories were compared with the numerical simulations 
on a preliminary basis in Ref. 2 and reasonably good agreement was noted. 

A more thorough review was undertaken and will be presented in the following 
pages. The difficulty in using the data is caused by the shortness of the 
histories presented and the rapid transients included in them. Although not 
of direct interest to the flow problems being considered, the inclusion of 
the transients provides a severe test for the code. 

The spatial dispersion of the steep gradients in the flow must be inter¬ 
preted to coincide with a point in space. The arrival of the shock wave at 
a point will be measured from the center of the dispersion and not from its 
leading edge. Similarly, if we are to compare space-time points, the pressure 
records must be adjusted to a point, rather than the average pressure they 
monitor. 

The paucity of information on how to approach the interpretation of the 
averaging of the applied pressure by transducers led to the inclusion of the 
analysis in Appendix B. This appendix is based on some referenced notes and 
should be useful to the reader in considering the averaging process in more 
detail. Using these techniques, a comparison of pressure monitored at the 
center of the gage versus the center of a comparably located cell could be 
undertaken. 

In the earlier comparison (Ref. 2), the gage-sensing element was 0.218 in. 
in diameter and the cells used were 0.125 in. square. As mentioned earlier, 
this geometric ratio gave identical "rise times. 11 This term is used to refer 
to the time it takes the sensed or calculated pressure (or other flow variables) 
to reach a step value in response to the incidence of a classical step shock 
wave, i.e., the time it took the shock front to cross the sensing element. 
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Although matched on a rise time basis, the pressure gage monitored an 
area 2.3 times greater than the cell of the numerical simulation. The impli¬ 
cations of the use of the averaging process involved for both simulations is 
described in Appendix B. Another consideration related to the use of experi¬ 
mental data for the evaluation of the numerous simulations is the measurement 
accuracy of the monitoring and recording systems. The systems had a 5 percent 
error in amplitude measurement and could respond in the range of 230 kHz 
without resonance (about a 3 ^sec rise-time capability in response to a step 
pulse). 

The present comparison is for the same problem, but 1/8- by 1/4-in. 
cells were used. The cells were positioned with their long sides parallel 
to the center line of the entry. The new cell size increases the resolution 
of the numerical simulations and represents about the coarsest grid employed 
in the numerical simulations for development of the conceptual model. The 
comparison of pressure histories for two gage positions is given in Fig. 6. 

Some deviations in the profiles can be noted in regions where sharp 
gradients occur. The experimental data were taken from reproductions in 
Ref. 4 and are estimated to be accurate to about 7 percent. Some of the 
resolution in the code was lost in using large printout intervals; however, 
the primary cause of the loss of the high-frequency components of the traces 
is the large cell size. 

In addition to losing some of the finer characteristics, there was also 
a phase shift (10 to 50 jxsec) in some regions of the histories. The phase 
shifts are due to the spatial dispersion of the wave in combination with the 
coarse grid sizes. Comparisons of the calculated data at selected other 
locations show equal or improved resolution. Since the original data were 
not available, further comparative work was not undertaken. The relatively 
good agreement for the coarse grid used provides sufficient confidence in 
the reliability of the code to reproduce the pressure distributions that 
correspond to the flow fields of interest. 
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A similar test comparison using the FLIC code was made for an entrance¬ 
way in Ref. 7. The results were similar — a good long-term pressure contour 
but a poor resolution of extremely short-term transients. The result is of 
interest to the shelter problem, because it supports the argument that the 
local pressure transients developed by shock waves moving through chambers 
with small orifices (tflO percent) are in effect really transients and do 
not substantially affect the ,f filling ,T process. If they were important, 
it would be difficult to explain why the long-term pressure contour of the 
numerical simulation compared favorably with the experimental data. 

A further review of the predictive accuracy of pressure histories deter¬ 
mined by codes as compared with test data from scale models will not be under¬ 
taken. Full-scale data will, however, be used in a later section. The 
remaining comments on the BRL model experimental data will relate to the 
smoke-grid experiments (Ref. 5). The geometry of the grid was given in 
Fig. 2. 

In the numerical simulation process, the shock wave’s discontinuity is 
spread over 3 to 4 cells ahead of the front. A "marked" particle placed in 
the flow will have a false path traced in Langrangian coordinates. The 
spreading of the discontinuity establishes a flow in the simulation ahead 
of the shock wave. The comparison of the path of marked particles in the 
simulation would not, therefore, correspond to the counterpart particles in 
the experimental flow. 

For example, Particle 1 in the simulation moved sooner and more rapidly, 
and at approximately 320 ^sec overtook the position of experimental Particle 2. 
From this point in time, Particle 1 of the experiment had the x,y,t path 
corresponding to simulated Particle 2. A method of comparing the experi¬ 
mental and numerical simulations under this set of conditions is to plot the 
velocity-time histories of each point. By establishing the x, y position of 
the particle as a function of time, the corresponding velocities of the 
experimental and numerical simulations can be compared. This procedure was 
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undertaken for the data presented in Fig. 7 and the results are plotted in 
Fig. 8. The plots provide velocity—time profiles for fixed points within the 
model chamber. The plotting method combined the spatial and velocity—time 
histories of each particle in the experiments and the simulation in order to 
determine the velocity—time distribution for fixed spatial points. 

This representation (Fig. 8) shows a better correlation between experi¬ 
mental and numerical data than the simple plots of the particle paths. The 
improved accuracy results from comparing velocities at like space-time points. 
The direct comparisons of particle trajectories included large leverages which 
magnified small errors. The magnification or leverage was caused by the strong 
spatiotemporal variations in the flow fields. Once a particle f s location was 
slightly in error relative to a fixed time, a very different history resulted 
because the flow-field gradients were so steep. 

The accuracy of the comparisons between experimental and numerical simu¬ 
lations ranged between 3 and 10 percent for most data points. This result is 
considerably smaller than the 10-to-15-percent error range estimated from the 
raw data in Ref. 2. The grid used in that simulation was fairly coarse (four 
cells were used to represent one-half the entrance opening) . Therefore, as 
a test of whether or not the accuracy of the simulation could be improved, 
the simulation was rerun using the grid shown in Fig. 9. Lagrangian marked 
particles were placed in the simulation, but the early motion of particles 
caused by the spreading of the shock wave again made comparisons difficult. 

To avoid the complex cross-plotting undertaken in the preceding work, 
the absolute velocities (the magnitude of the velocity vector) of two particles 
in the experiments were compared to those developed by the numerical simula¬ 
tions. The spatial position of each particle at a given time was determined 
from the experimental data (Ref. 4) . An interpolation for the absolute ve¬ 
locity at an equal position in space-time was made for the numerical simula¬ 
tion. The results are presented for two particles in Fig, 10. 
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Fig. 8. Velocity Time Profiles for Fixed Positions in the BRL Model 
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Fig. 10. Velocity—Time Profiles 
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An increase in the deviation between numerical and experimental data 
is evident. The time grid did pose some problems in the interpolations in 
location and time. A 30-^sec variation in the profile is plotted in the 
upper portion of Fig. 10, and a large change in the deviations between the 
experimental and numerical curves can be observed. Linear interpolations 
were used to determine the axial and radial velocities at a space-time point 
in the numerical simulations from which the absolute velocity was calculated. 

A more sophisticated and accurate interpolation procedure would have been 
preferable, but the time involved in establishing the procedure and intro¬ 
ducing the data could not be justified under the limited scope guidelines of 
the project. 

The presentation of the techniques for interpreting pressure—time his¬ 
tories outlined in Appendix B discussed some of the problems associated with 
determining the values of variables at specific space-time points. A direct 
quantitative application cannot be made because the assumptions made regarding 
the relationship of wave shape and cell size do not hold for the general cases 
of the simulations. The simulations, in addition, present data for specific 
points and not areas, as was the basis for the methods presented. The inter¬ 
polation could use the general guidelines established; but, more importantly, 
the material did provide some insight into the problems of interpreting 
spatially and temporally varying data. 

A simpler procedure than plotting the velocity-time histories was used 
for further comparisons of the data. The vector values of velocity (magnitude 
and direction) were established from the experimental and numerical data for 
the experimental space points at three times. Time was measured from inci¬ 
dence of the shock wave at the inner surface of the entrance. Time interpo¬ 
lation was avoided by selection of the points of time comparisons, and a 
variation of under 7 /jt sec between comparisons was achieved. Linear spatial 
interpolation was still used to provide correlations between corresponding 
points. Linear interpolations were also used by Coulter to develop the 
experimental velocity vectors from the particle trajectory data. 
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Tables 1, 2, and 3 present the data for the 20 grid points shown in 
Fig. 2. The experimental values used in the evaluations were the corrected 
values from the errata for Ref. 4. In Table 1 (213 jLisec) the deviations 
between numerical and experimental data range from 1 to 30 percent. The nu¬ 
merical simulation showed higher velocities along the centerline of the door¬ 
way and lower velocities away from the entrance. The vectors had a deviation 
range of 0 to 18 deg, with 85 percent being within 8 deg of each other and 
95 percent within 11 deg. 

At 450 ^sec (Table 2) , the variation in the values of vector angles in¬ 
creased slightly, as did the variations in flow magnitudes. The trend con¬ 
tinued and increased at 750 ^sec (Table 3). For the full time range, the 
centerline and boundary flow vector angles agreed fairly well. The centerline 
numerically simulated flows remained above the experimental values and the 
off-axis flows continued below the experimental values. 

The first inclination is to fault the numerical simulations; however, a 
closer examination of the experimental data shows some inconsistencies. Point 
5 in Table 2 has an unrealistic magnitude, and the variations in vector angle 
for points 9 and 10 are inconsistent with the flow field geometry. In Table 3, 
the flow magnitudes at experimental points 6, 19 and 20 appear unrealistically 
large. The vectors at experimental points 10, 14 and 20 are inconsistent with 
the flow-field geometry. These inconsistencies, coupled with the errors 
caused by small temporal and spatial variations, prompted a brief investiga¬ 
tion of the accuracy of the experimental data. 


The first area investigated was the data collection and reduction process. 
The camera framed a picture of 1-^sec duration every 38 jusec. The zero time 
had to be computed as the shock wave passed the reference plane between frames. 
Due to the relative invariance of shock front velocity (1282 fps incident and 
1132 fps at sound speed) an error of only a few microseconds should have re¬ 
sulted. More error would be attributed to the reading process, since grid 
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Magnitude of the velocity vector of angle 6 and x component u and y component 





FLOW VECTOR COMPARISON AT TIME = 406 jusec 
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Corresponds to smoke grid intersection as defined in Fig. 2. 

Magnitude of the velocity vector of angle 9 and x component u and y component 
Sign conventions shown in Fig. 2. 




FLOW VECTOR COMPARISON AT TIME = 713 /j.sec 
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Corresponds to smoke grid intersections as defined in Fig. 2. 

Magnitude of the velocity vector of angle 0 and x component u and y component 
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intersections must be read. The intersections are formed by finite-width grids 
which distend in time and space. The grid intersections are small blocks simi¬ 
lar to the blocks formed by the grids from which the densities were calculated. 

The edges become more diffuse with time, hence each intersection point 
may not represent the same particle over a period of time. The extensiveness 
of the turbulent flows and the boundary layer growth along the viewing windows 
in the later portions of the flow field would add to the diffusion process. 

The local, as well as the regional, momentum transfers will also affect the 
experimental grid locations. 

The averaging interpolation process in evolving the numerically simulated 
data was found to be a highly sensitive factor. The linear method chosen did 
not always present the most accurate result, but it did provide a consistent, 
simple and unbiased value. It is expected that the averaging process for 
evolving the experimental velocities and vectors would have encountered similar 
difficulties. An error analysis of the averaging processes would be a fruit¬ 
ful exercise if further evaluations are undertaken. 

Because of constraints on the scope of the project, a less exacting and 
simpler method was chosen. The front and rear smoke grids experiments for the 
same model at equal incident conditions were collated as shown in Figs. 11, 12, 
and 13. The data at equal times should be consistent for both experiments (a 
small time bias of 18 fj ,sec exists in the presented data and the timing differ¬ 
ence between the data presented in Ref. 2) . The velocity magnitudes and vectors 
are presented as a function of location in the flow field. 

In Fig. 11, the region in the x direction (labeled A) shows an inversion 
of velocities at the interface between the front and rear model data. Some 
small variations between the vector directions and the normal to isovelocity 
contours can be observed. A variation of 20 to 30 fps at 150 fps and a vector 
variation of 10 to 20 deg appear to be reasonable estimates for average values. 
The flows in Region B are greater than those expected from other flow data. 

For example, the 257-fps value exceeds the flow velocity of the incident shock 
wave and is off axis from the axially directed inflow. 
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Fig. 12. Experimental Flow Contours at Time = 405 /nsec 









UR5 


755-3 


In Fig. 12, Regions A, B, C, and D show some unusual variations in 
either flow magnitude, direction, or both. The data show a slightly smaller 
variation between front and rear flow magnitudes (relative to Fig. 11) , but 
indicate equal or slightly increased angular variations. Region A shows a 
wide range of variation in flow velocity magnitude and, to a lesser degree, 
direction. Region B shows some rather unusual changes in flow magnitude. The 
flow magnitudes in Region C appear considerably higher than the surrounding 
regions. In Region D the flow vectors appear to have some inconsistencies 
in direction. 

Figure 13 shows an even greater increase in the variation between front 
and rear data of the BRL experiments. Region A in Fig. 13 indicates flows 
inconsistent both in magnitude and in direction. The flow vectors shown in 
Region B have reversed flow magnitudes compared to surrounding regions. Re¬ 
gion C contains a flow vector with a magnitude and direction which is inconsis¬ 
tent with the flow field. Region D includes a wide range of variations, which 
may indicate a pronounced turbulence or the limits of experimental accuracy 
at this time in the flow history. Region E indicates a gradient in the flow 
magnitude that appears inconsistent with the flow field at this point in its 
history. The range of variations appears to be 15 to 35 percent in magnitude 
and 10 to 40 deg in direction. 

A comprehensive statistical analysis of the data was not within the pro¬ 
ject scope; hence a simpler approach was taken to provide some quantitative 
estimate of the variations. The method selected plots the angular deviations 

(A 9) and the ratio of vector velocities (R ) , numerical to experimental, 

ne 

presented in Tables 1, 2, and 3 on normal and log-normal probability graph 
paper. Assuming the errors occur randomly about some mean value, the data 
points should form straight-line segments. Each segment would represent an 
error distribution in a region of the flow. 
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The data plotted from Table 1 indicated three distinct distributions for 

the ratio of the magnitudes of velocity vectors (R ): one for the inflow 

ne 

region (points 1, 2, 3, 4, and 8), and a limited one for points 13, 6, 7 
and the common distribution for the remainder of the points. The mean for the 
velocity ratio first set was 1.28; for the second set (points 6, 7, 13) the 
mean was 1.02; and for the third set the mean was 0.82, The vector direction 
variations (A0) showed two regions, one containing points 1, 2, 3, 4, 5, 6, 10, 

11, 17, 18, 19, 20 and the other the remainder of the points. For the flow 
direction deviations the means were 3 deg and 8 deg, respectively. 

The least variation in flow direction between numerical and experimental 

data was in regions where the flow was the strongest, or where the flow was 

along a wall. The values of the R (numerical-to-experimental flow-magnitude 

ne 

ratios) showed two distinct flow areas, with an increase in variation between 
the inflow region and the bounding wall. To help understand why the error 
increased, iso-R^, iso-AO, and isovelocity contours were plotted. 

The plots reduced the indicated variation between numerical and experi¬ 
mental simulations in the central region of the off-axis flow. If corrections 
of the flow variables were developed from the R contour values, the average 
velocity ratio variation rose to about 0.93 for the off-axis region. The in¬ 
flow ratio was not altered. The distributions in both zones appeared to be 
random. 

The same procedure was repeated for the data in Table 2. A similar two- 
segment distribution for the angular deviations was acquired. The mean of 
the variations increased for the distributions from 2 deg to 4 deg and from 
8 deg to 20 deg. 

The velocity ratio distribution changed slightly. The inflow zone mean 
increased to 1.33, as did the variance. The second zone consisted of points 

12, 7, 13, 9, with points 8, 3, 11 forming its lower bound. The variance in 
this zone increased, and the mean was 0.92. The third segment of the plotted 
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velocity data contained points 6, 10, 11, 14, 16, 17, with a mean of 0.55, 
and the fourth segment contained points 15, 18, 19, 20, with a mean of 0.30. 

Each set of points was roughly related to zones of flow in the chamber, but 
a strongly decreased consistency in the relationships was noted. 

If adjustments of purported inconsistencies in Fig. 12 were made, a 
complete alteration in the above results was noted. Mean values of the R 
ratios began to approach 0.80 to 0.85, but a change in the values along the 
inflow axis was again absent. A gradient did occur between the inflow value 
and the flow parallel to the entrance axis that corresponded to the turbulent 
region of the jet. 

The data from Table 3 shows a continued decrease in correlation between 
experimental and numerical simulations. The distributions of velocity ratios 
increased, as did the variations in flow direction. The mean values were all 
higher, with the value of the inflow ratio R n0 approaching 1.90 at particle 
position 4. A marked increase in variance occurred in this region. 

The distribution of the data plotted from Table 3 was completely random 
in character. Particle position did not seem to be correlated with the re¬ 
sultant distributions. The isovelocity and isoratio contours did not provide 
additional insight into the nature of the distributions. The distributions 
formed into approximately the same two zones as did the data from Table 1, but 
with means of 1.65 and 0.75. If adjustments for the inconsistencies noted for 
Fig. 13 are made, the mean of 0.75 cited above for the off axis flows increases 
to a value in the range of 0.85 to 0.95. 

The following conclusions were drawn from investigations. The experimental 
flow measurements increased in error with time, with strong gradients in flow 
magnitude and direction, and with decreased velocity. Errors of 20 percent in 
magnitude and 20 deg in direction are not uncommon and can be significantly 
greater at late periods in the flow history. An exhaustive error analysis is 
needed to fully interpret the data. 
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Relative to the numerical simulations, some observations on their accuracy 
can be made. The code does introduce some error in the turbulent inflows of 
the models tested. These turbulent regions have velocities 5 to 15 percent 
less than predicted by the numerical simulations. The overall flow patterns, 
both long and short term, are similar in character and magnitude for both 
simulation techniques. 

The code produced higher inflows in the second simulation discussed. 

This deviation is attributed to the difference in entrance conditions between 
the two numerical simulations and the experiment. The earlier simulation, 
in which excellent agreement of centerline velocity profiles was reached, 
more closely resembled the upper portion of the BRL model; the second simula¬ 
tion more closely resembled the lower half (see Fig. 14). These variations 
in the two simulations were the result of using axial symmetry to save on com¬ 
puter time and memory and to increase resolution. Upstream pressures were 
not available from the model experiments, but a large-scale test (see Appendix 
C) similar to the first simulation showed a decrease in reservoir pressure of 
3 psi. 

The magnitude of the flow variation caused by the difference in reservoir 
pressure would constitute a portion of the variation noted. The increases 
in pressures on large reflecting surfaces and reservoirs of larger volume were 
noted and discussed for chambers in Ref. 42. This effect could be responsible 
for the higher inflows noted in the second numerical simulation, but it does 
not fully explain the higher velocity flows along the entrance axis later in 
the flow cycle. 

Viscous and turbulent losses during the inflow and the filling cycle are 
to be expected, particularly in the smaller scale models. The comparison of 
numerical and experimental data supports this argument. The flow magnitudes 
of the simulations were about 10 percent higher along the inflow axis and 5 
to 10 percent lower in the remainder of the chamber. 
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Although these estimates of numerical simulation accuracies are not as 
precise as desired, further evaluations do not seem warranted. There is diffi¬ 
culty in providing precise comparisons between the parameters monitored in 
small-scale models and the numerical simulations. Small errors become signifi¬ 
cant. It is difficult to ascertain the absolute accuracy of experimental data. 

The intent of the original work plan had been to secure the experimental 
digitized data from the original cards or tape used to develop the appendices 
of Refs. 4 and 5. This plan proved infeasible; hence the experimental data 
used in the comparisons either had to be reproduced in a format compatible 
with the Office of Civil Defense's CDC 3604 computer or had to be worked with 
manually. Either approach was severely limiting to the evaluation and compar- 

j|c 

ison procedure. 

For these reasons, further evaluations of the code using experimental 
data from scaled experiments were not undertaken. The shelter program's inter¬ 
est in large-scale data also suggested that more emphasis be placed on evaluation 
of the large-scale situation. Two sources of data were used: field test data 
and URS shock tunnel tests. The full-scale field test data used for evaluation 
were those presented by G. Coulter in Ref, 5. Though not precisely two-dimension¬ 
al, the data were considered to adequately approximate the two-dimensional case 
in sufficient detail to test the numerical simulation. 


* Investigations under this study indicated magnetic tape conversion or the 
development of punch cards from published tables to be too costly an under¬ 
taking for studies such as this work unit. Future OCD requirements would 
be best met if all new digitized experimental data is in a format compat¬ 
ible with the OCD computer. 
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COMPARISON OF NUMERICAL SIMULATION WITH FULL SCALE EXPERIMENTAL RESULTS 

The chamber geometry for the large-scale test is presented in Fig. 15 
and the geometry of the numerical simulation in Fig. 16. The short clearing 
time of the structure and the entrance near the side of the building (and 
away from the stagnation region in the center of the front wall) reduced 
effects of the three-dimensional flows. In the two-dimensional case, the 
reservoir pressure would be slightly greater due to the absence of the upward 
expansion of flow over the roof of the structure. The higher reservoir pres¬ 
sure would result in higher inflow velocities and a higher internal pressure. 

The monitored free-field pressure was used as an input to the numerically 
modeled situation. The input took the form of a calculated pressure-distance¬ 
time profile whose propagation over the measurement point which would correspond 
to the experimentally monitored pressure—time history. The results of the 
simulation are plotted in Fig. 15, as were the experimental data as taken from 
Ref. 5. To save computer time, the numerical simulations were only carried 
out to 60 msec, which is well into the outflow stage. No significant alter¬ 
ation in accuracy between the numerical-to-experimental profiles would be 
expected after this point. 

The numerical simulations provided excellent agreement (better than 
20 percent) for the first 20 to 25 msec. At this time, the three-dimensional 
effects began to be sensed at the internal gage locations, and the error in¬ 
creased 10 to 20 percent. An alternate view of the error is to consider it to 
be a 5- to 10-msec lag in the decay of the pressure. This agreement is better 
than that with the small-scale models discussed earlier and is well within 
the needs of this study, which is primarily interested in the more nearly 
two-dimensional basement shelters. Based on this comparison, the accuracy of 
the code in predicting pressure distributions in shelters would appear to be 
well above 90 percent. 
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Fig. 15. Comparison of Numerical and Experimental Full Scale 
Pressure Histories 
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Fig. 16. A Two Dimensional Simulation of a Free Field Test Case 
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An evaluation of full-scale test is presented in Appendix C. Due 
to the lack of other good full-scale data, it was decided to make use of 
the data in this study. The URS version of the FLIC code was used to 
numerically simulate the tests. Because of the length of the presenta¬ 
tions and comparison of the results, the material was placed in Appendix 
C. The comparisons of test results support the conclusion that full-scale 
numerical and experimental simulations seem to agree to within a 90-per¬ 
cent absolute accuracy on both overpressure and dynamic pressure predic¬ 
tions , 

The main difficulty encountered in comparing the numerical and experimental 
simulations was the manipulation of the available experimental data into a use¬ 
ful format. The accuracy of the distributions of two parameters were of prime 
interest: pressure for structural loading and dynamic pressure for human 

injury considerations. Pressure correlations are direct measurements and, if 
a common time base for comparisons can be established, comparisons are simple. 

As was noted earlier, the small model data did provide some difficulties in 
matching time scales. The problem did not reappear in the full-scale experi- 
mental-to-numerical comparisons. It was found that, in general, the experimen¬ 
talists did not provide sufficient error analyses with their pressure data to 
allow precise comparisons to be made. As the variation between experiment and 
numerical simulation dropped below 10 percent, it was difficult to fault the 
numerical simulations, inasmuch as exact pressure-amplitude and time-scale 
error values and their distributions were not available. 

The flow measurements underscored the problem more clearly. The measure¬ 
ments in the smoke grid experiments were presented to the thousandth of an 
inch (Ref. 5) and the time scale to the 100 nanoseconds. Even if this precision 

really did occur, two views of the same process did not overlay accurately. The 

■ * 

preliminary error analysis accompanying the URS preliminary shock tunnel test 
data showed where large errors could be produced from basically accurate data. 
The question repeatedly arose as to how accurate were the data against which 
the numerical simulations were being compared. 
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Several factors noted in these comparisons will be presented here; other 
factors will be covered in later portions of the report. Although the exper¬ 
imental model data showed a pronounced spreading of the inflow through tur¬ 
bulence, this was not observed in the large-scale tests. The scale-model tests 
also showed an increase in turbulence with decrease in entrance size. The 
tests of the 1/8-in. entry showed an interference between opposite sides of 
the inflow pattern and the establishment of a Von Karman vortex street. These 
combined factors indicate that perhaps a scale effect may exist. 

A marked effect on the inflow is exhibited by the reservoir conditions. 

As the entry cross section becomes of appreciable size (10 to 15 percent) in 
relation to the reservoir cross section, an appreciable change in reservoir 
pressure occurs. For a 19-percent opening in the URS shock tunnel tests, 
approximately a 29-percent decrease in reservoir pressure was monitored. This 
change in reservoir pressure was more rapid than the time normally ascribed 
to clearing phenomena. 

Questions as to the distribution of kinetic and potential energy in the 
inflow were raised. Significant flows were monitored in the chambers after 
the filling process (as evidenced by the pressure increase) had terminated. 

The momentum generated by the inflow process was sufficient to maintain the 
flow, and form the circulation referred to in preceding work units. The 
review indicated the phenomenon was strongly associated with the redistribu¬ 
tion of energy, and it appeared to be a critical factor in understanding the 
flow process. 

The evaluations could not establish any significant evidence that the 
numerical simulation of full-scale situations would have errors over 10 per¬ 
cent in either overpressure or dynamic pressure. The evaluations did seem to 
indicate that if any errors did occur in the simulations, they tended to 
produce greater values of both variables. As to the predictive validity of 
the numerical simulations, some doubt remains regarding their possible over¬ 
emphasis of the long-term flow in the chambers. This overemphasis would be 
expected in the absence of turbulent or viscous dissipation in the numerical 
simulations. 
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The remaining inconsistency of concern to this evaluation was the higher 
inflow velocities developed by the numerical simulations in the central and 
aft regions of the chamber particularly where the inflow turned and moved 
along the rear wall. The variation could be due to Galilean invariance of 
the numerical techniques, as discussed earlier; however, the pressure histories 
monitored on the rear walls of the full-scale tests (field and URS shock tunnel) 
did not indicate similar variation. The gages on these tests were exposed to 
stagnation pressures for the turning flow. The indication is that if an error 
does occur in this region, it is highly localized. The experimental data from 
the scaled models were slightly inconsistent at this point in time. Therefore 
it is hard to pinpoint where the error lies. With the agreement of the simu¬ 
lation with large-scale data, there appears to be no reason not to use the 
results of the simulations in these regions. 

The investigation undertaken could not establish a reason for not using 
the two-dimensional hydrodynamic code in the OCD shelter research program. As 
new data are developed, this conclusion should be reviewed, particularly with 
respect to long-term flows. The predictive accuracy of the numerical simula¬ 
tions with respect to pressure distributions were particularly good. No evi¬ 
dence could be developed that indicated any significant variations between 
code and experiment relative to the information on structural loading required 
by the OCD shelter research program. The overall usefulness of the codes to 
the shelter research program appears to be more a function of the factors of 
utilization, cost, and data display. 

UTILIZATION, COST, AND DISPLAY OF NUMERICAL SIMULATIONS 

The term "numerical simulation" has been used to discuss the code outputs. 
The choice was selected to indicate the close resemblance, in terms of the data 
acquired, between using the codes and undertaking experimental studies. The 
process of planning and executing a simulation follows a procedure identical to 
establishing and executing an experimental test plan. The main exception is 
that the data reduction process entailed in experimental work is absent 
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Under this assumption, it would be useful to compare the features of both 
approaches to simulating the shelter environment during a nuclear attack. One 
of the main features favoring the experimental approach is credibility. Where 
numerical and analytical studies suffer from a credibility gap, the experimental 
results tend to be unquestioned. In the preceding review some difficulties 
were noted in both techniques. These difficulties will be reviewed to form a 
background for the discussion of the usefulness of the numerical simulation 
techniques. 

One of the problems noted in the preceding discussion was the absence of 
error analysis on, and verification of, experimental results in the flow studies 
as was underscored by the comparisons against the numerical simulations and 
their unbiased consistency. The numerical simulations covered a broad spectrum 
of geometric scale, which underscored the fact that the shelter research pro¬ 
gram is interested in what happens in the full-scale prototype situation and 
not in the simulations. The simulations, both numerical and experimental, are 
the means to the end, not the end itself. 

Once stated, the above is so obvious that it appears too elementary to 
merit consideration. However, reviewing the available material, most studies 
did not appear to have a strong problem-oriented focus. This situation appears 
to have been caused by having to use small-scale experiments to study the 
physical process due to the expense and large yields resulting from full-scale 
testing. An example of the potential conflict can be found in the chamber¬ 
filling process. The chamber-filling process, based on the volume-to-area 
ratio model of Ref. 42, is nonlinear; yet the validity of linear scaling of 
the flow processes has not been questioned. 

The difficulty of acquiring valid and comprehensive flow data must also 
be reviewed. The experimentalist selects monitoring locations which he feels 
will provide a true picture of the phenomena being investigated. Constraints 
governed by cost, available instrumentation, and size (in scaled experiments) 
limit his choice in varying degrees. The combination of the above factors can 
produce some misleading results. 
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An example of these misleading results can be found in Ref. 43, in which 
smoke-paricle experiments were used to visualize flow in chambers. These ex¬ 
periments observed the motion of a single smoke grid initially positioned 
across the door entrance of the models shown in Fig. 17. The trajectory of 
the grid was used to establish the "maximum flow" in the chamber. These ex¬ 
periments are similar to those in Refs. 4 and 5, which used a grid of smoke 
particles (Fig. 2) to establish a more general flow pattern. 

The maximum velocities reported are not consistent with either the data 
presented in this report or in Ref. 2, or from considering the mass influx 
necessary to fill the chamber in the correct "fill time." The velocities pre¬ 
sented in Ref. 43 are considerably lower than those indicated in the other 
sources. The variation is caused by the time required for the inflow to ac¬ 
celerate to its maximum velocity. During this period, the smoke particles 
move away from the high-flow zone. The conclusion in Ref. 43 (on the spatial 
flow distribution) is drawn from observations of the particles traversing the 
inflow zone, but does not account for the strong temporal variations which 
exist in the flow field. 

Both sets of experimental work reviewed above were extremely well done and 
produced good and consistent information within the experimental limitations. 
These limitations precluded long-term observations of the spatial distribution 
of the flow field. This limitation, in turn, did not provide the needed over¬ 
view to develop an entire review of the flow fields generated. These comments 
are not presented as an argument against experimentation, but rather as re¬ 
minder that the data may not be complete or infallible, particularly in cases 
where strong temporal and spatial variations of multiple parameters exist. 

Utilization 

The numerical simulation technique does provide an advantage which cir¬ 
cumvents the problem of having to define monitoring locations. The code must 
calculate the complete set of parameters needed to define the flow field at 
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DIRECTION OF INCIDENT SHOCK WAVES 




SHOCK TUBE WALL 



SHELTER WITH UNPROTECTED DOOR 




SHOCK TUBE WALL 


SHELTER WITH PROTECTED DOOR 
(Reproduced from Ref. 43) 


Fig. 17. Orientation of Shelter Models with Respect to the Incident 
Shock Waves for Tests in Ref. 43 
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every cell. The complete set of data for all cells in any form desired can 
be provided in the output. The cost of the procedure is dependent on the 
exact output procedure, whether output is taken from all the cells or just 
some portion of them, and on the particular computer involved; but in general 
the cost range is about 5 to 30 percent of the total program running time. 

The second advantage of numerical simulation is that the output is in 
the form of final data, not raw data which must be reduced to a useful form. 

The absence of the data reduction step has extremely attractive economic aspects. 
By placing the output of the code onto magnetic tape, data searches, data 
ordering, and data plotting can be done quite simply and economically. The 
overall effect is a flexibility, thoroughness, and rapidity which experimental 
techniques cannot often approach. 

Two basic arguments can be presented against numerical simulations. 

These concern the credibility of results, since the simulation is not "physi¬ 
cally real,” and consideration of cost. The credibility objection can often 
be valid, but the code can be checked against limited experimental data for 
verification. The difficulties of a complete verification of numerical with 
experimental simulation were discussed earlier in considering code accuracy. 

Part of the difficulty is that the experimental data relating to flow param¬ 
eters is not complete. However, since data are now available for study that 
were not obtainable from preceding studies, some estimate of the relationship 
of numerical simulation to experimental data can be inferred. 

Another aspect of credibility must also be probed, i.e., the fact that 
experiment, like the code, is still a simulation of the prototype (full-scale) 
situation. The situation could exist in which the small-scale experiment could 
not be simulated with the accuracy required for a simulation of the prototype 
situation. The comparison of code to experiment made earlier in this section 
did indicate the code was slightly more accurate in the large-scale case than 
in the modeled cases (although in both cases the accuracy was sufficient for 
the purposes of the shelter research program). Based on this reasoning, the 
credibility argument against numerical simulation, as applied to this OCD 
shelter research program, is weak. 
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Cost 


The second argument often used against the utilization of numerical sim¬ 
ulation concerns cost. The only objective evaluation of cost available is to 
compare it to the other form of simulation — experiment. A comparison of the 
procedures involved in acquiring data through each of the simulation procedures 
is presented in Table 4. The assumption is made that the experimental facility, 
the data acquisition system, the experimental data reduction system (including 
the computer programs for data reduction), the computer, and the numerical 
code are available for project use on a per-hour basis. The table indicates 
that the first and final steps are similar for both simulation procedures, 
and it will be assumed that they represent equal costs. 

Under these rules the comparison lies in the intermediate stages. These 
stages (as outlined in Table 4) are formed by experimental steps 2, 3, and 4, 
which correspond to numerical step 2. From past experience, step 2 in either 
simulation procedure costs approximately the same on either project. The 
experimental costs would be greater than the costs of the numerical simulation 
by an amount equal to the cost of experimental steps 3 and 4. The exact dollar 
cost would depend upon experimental conditions, such as size, data acquisition 
needs, number of times the model will be tested, and the complexity of the 
experimental design. 

The next stage in the comparison is the actual performance of either 
simulation. Each of the shelter chambers numerically simulated for OCD by 
URS have cost in the range of $30.00 to $120.00 at standard commercial rates. 

If the pretest cost of the experimental facility and data acquisition system 
is neglected, the per-test cost is the project charge of the experimentalist 
and his test crew. On the basis of standard overhead and salary rates, the 
project cost of a senior experimentalist would be $20/hr and a senior technician 
$12/hr for a combined rate of $32/hr. Depending on the size of the facility, 
more technicians might be used, thus raising the per-hour cost further. The 
cost per experiment for either size test group would be comparable or greater 
than the cost of numerical simulation without the overhead costs of the ex¬ 
perimental facility being involved. 
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Table 4 

AN OUTLINE OF NUMERICAL AND EXPERIMENTAL SIMUIATION PROCEDURES 

Experiment Numerical Simulation 

1. Design the experiment 1 # Design the experiment 

2. Design the data model 2. Develop the input data 

3. Fabricate the model 3. Perform the simulation 


4. Design and set up the data 
acquisition system 

5. Perform the experiment 

6. Reduce the data 

7 # Perform an error analysis 
of the experiment 

8. Data presentation 


4. Review the output for errors or 
inconsistencies 

5, Data presentation 
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Another pertinent comparison would be the use of the OCD Olney Computer 
Facility for the numerical simulation. This facility is provided to projects 
on a time-available basis. The cost to the project of each simulation, in¬ 
cluding sending a person from the most distant location in the country, ranges 
between $5 and $20 for batch runs, depending on the number of simulations in 
each batch. The cost of the performance stage of the simulation will, in either 
situation, be equal to or less than the cost of experiment. 

The last items listed in Table 4 compare experimental steps 6 and 7 with 
numerical simulation step 4. The data reduction process is not necessary in the 
numerical simulations; thus, the added cost of the error analyses and consist¬ 
ency checks are equal in both simulations. The computer costs of adequate data 
reduction is similar to the computer costs of the numerical simulation; hence 
a significant difference would exist between simulations. If all these cost 
factors are summed, numerical simulation has an intrinsically favorable cost- 
benefit ratio. If we further define the ratio to account for the amount of 
information per simulation, the more complete data from the numerical simulations 
still further increases the cost-benefit ratio in favor of numerical simulation. 

Another, more powerful, approach would combine numerical and experimental 
simulation. In this case the greater portion of the simulation would be per¬ 
formed numerically, reinforced by experiments to substantiate the numerical 
results. The combination provides the credibility of experiment with the 
thoroughness and cost-saving features of the numerical simulation. 

Display 

The combined approach still contains the data display difficulty found 
in both simulation techniques. Both simulation techniques simply provide data, 
not mathematical and physical understanding of the flows, which reduces both 
simulation techniques to the status of tools for the researcher, i.e., they 
are means to an end rather than end product. The inability of many studies on 
shelter flows to display and use available data stresses the lack of under¬ 
standing of the basic processes involved. The code, like experiment, provides 
data for theoretical or empirical analysis, which is the basis of the shelter 
research program. 
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The problem with data display is primarily due to the wide variety of 

shelter geometries and input conditions possible, to the strong spatio-temporal 

✓ 

variations, and the many flow variables needed to describe the chamber flows* 

In the chamber-filling studies (Refs. 2 through 5, 42, and 44), the methods 
for predicting a gross pressure-time history in a chamber became complex and 
the authors finally resorted to a series of simple computer codes. At present 
this is still the major shortcoming of both simulation techniques. 

The cost-benefit ratios of the numerical simulations were further improved. 
The computer codes allow flow out of three of the four sides of the computa¬ 
tional mesh. The fourth side, the input flow boundary, can only have inflow 
and cannot accept the incidence of a change in the flow moving upstream. The 
reflected shock wave from the chamber in the simulation, upon reaching the 
input boundary, terminates the calculation period. To provide acceptable test 
durations, a long distance between the chamber and the input boundary must 
exist. 

If equal cell sizes are used, then large computational meshes must be 
used. To reduce this requirement, variable cell lengths can be instituted. The 
cells become increasingly elongated in moving from the flow shelter interaction 
area toward the input boundary. Since the upstream region of the grid is not 
of interest, the coarseness of the grid is not of importance. The use of this 
technique effectively minimizes computer memory requirements. 

The use of variable cell size can be extended further. In regions of 
steep gradients and intense flow, the need for high resolution exists; but, if 
an adequate grid representation is used, computer memory requirements become 
excessive. To circumvent this difficulty, a doubly variable cell geometry 
can be used. With this approach, a 3000-cell mesh can provide quite high 
local resolution and still avoid long computation times. A 3000-cell grid 
is compatible with the 53-K bit storage available on the OCD Olney computer. 

The difference between the effect of a doubly variable grid size and a 
simulation with a fixed grid must be negligible if the technique is to be 
useful. The profile in Fig. 5 indicates the transmission of a shock wave is 
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not significantly perturbed by the variation in grid size; hence a test was 
run to show that long-term effects were similary unaffected. An existing 
test of constant grid size run on the OCD computer in late 1968 was selected 
and rezoned as shown in Fig. 18. The rezoned shelter was run on a GDC 6600. 

The results of both tests are presented in Table 5. No significant variation 
was noted, and the data provide a good indication of the repeatability of re¬ 
sults. It should be noted the runs were on different machines and at different 
facilities (commercial and government) , which also indicates the small vari¬ 
ability of running on different equipment. 

The main drawback of the above approach is that the computational time 
step is a direct function of cell size. The smallest zones then tend to con¬ 
trol the time step for the entire grid, which results in long running times. 

The code could be rewritten to allow calculations in small time steps. At 
certain integer steps, the next set of largg zones would be included. Then, 
after a series of these computational sequences, the next larger set of zones 
included, and so on. For the immediate needs of the shelter research program, 
this approach is an unattractive choice because of the potential reprogramming, 
debugging, and proofing costs. 

Still another procedure might be considered for future use —that of re¬ 
zoning. The flow pattern in later flows in the filling of the shelter by the 
air blast environment does not require the high resolution of the earlier 
flows. The grid could be rezoned to coarser zones to take advantage of this 
feature of the flow. The rezoning of the computational grid would reduce the 
number of cells and would also increase the computational time step, since 
the cells would be larger and the flows lower, all of which would reduce the 
running time. After reviewing this procedure, it was also discarded since no 
significant addition to the project f s objectives would have resulted. 

The last item to be considered under the computer code evaluation is 
the effect of the grid size on the predictive accuracy of the numerical simu¬ 
lation. Of importance to the shelter research program is the cost saving 
coarser grids would provide in reducing computer running time. A possible 
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Table 5 

COMPARISON OF MULTIPLY VARIABLE CELLS 
WITH NON VARIABLE CELL SIZE 




Axial Velocity 

Radial Velocity 

Overpressure 


Time 

(fps) 

(fps) 

(Psi) 

Point 








(msec) 

5 cell 

10 cell 

5 cell 

10 cell 

5 cell 

10 cell 


n 

459 

464 

- 10.6 

-11.5 

4.86 

4.95 



613 

610 

16.6 

13.3 

3.08 

3.11 

i 

H 

613 

608 

16.7 

17.8 

3.92 

3.91 


40 

430 

432 

10.5 

13.2 

7.52 

7.45 


10 

480 

484 

-41.8 

-46.3 

3.33 

3.42 

o 

20 

538 

541 

-85.0 

-83.5 

1.53 

1.68 


30 

501 

503 

-83.0 

-83.9 

2.54 

2.46 


40 

233 

239 

- 39.4 

-44.7 

7.10 

7.03 


10 

242 

243 

14.8 

14.4 

2.56 

2.53 


20 

517 

514 

21.3 

22.5 

2.67 

2.61 

J 

30 

583 

574 

16.2 

16.0 

3.92 

4.03 


40 

609 

602 


3.40 

6.77 

6.82 


10 

231 

232 

74.2 

73.4 

2.55 

2.54 

A 

20 

449 

448 

84.8 

87.3 

2.59 

2.61 

4 

30 

453 

449 

66.0 

66.8 

3.89 

4.00 


40 

419 

418 

5.63 

4.08 

6.45 

6.54 


10 

60,3 

62.6 

3.42 

3.53 

1.01 

1.05 

r 

20 

173. 

176. 

29.2 

3.18 

3.49 

3.64 

o 

30 

129. 

128. 

11.7 

12.0 

7.57 

7.56 


40 

294. 

292, 

13.0 

11.5 

7.79 

7.88 



51. 

52.9 

9.8 

10.1 

.89 

.93 

c 

20 

171. 

170. 

17.48 

19.6 

3.76 

3.77 

o 


100. 

102. 

75.7 

74.1 

7.41 

7.42 


40 

219. 

223. 

81.6 

71.6 

7.84 

7.83 


10 

-.98 

1.0 

-.006 

-.0065 

0 

0 

17 

20 

9.8 

10.4 

5.4 

4.85 

3.83 

3.96 

7 

30 

21.1 

22.3 

8.78 

7.0 

7.65 

7.60 


40 

43.2 

45.1 

15.1 

10.6 

8.98 

8.85 


10 

-.69 

-.71 

-.188 

-.153 

0 

0 

8 

20 

9.4 

10.3 

31.35 

30.6 

3.87 

4.02 

30 

18.7 

19.0 

52.8 

53.4 

7.53 

7.58 


40 

31.4 

32.3 

45.6 

66.5 

8.22 

8.34 
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additional feature is that they may provide some insight into the feasibility 
of writing a small coarse-grid code to provide the general shelter research 
program and the practicing architectural engineer with a way of establishing 
load and flow information. This approach is attractive in that it limits the 
data display problem, provides for a wide range in geometric and input variation, 
and offers potential simplicity for the user. The key to its feasibility is 
the coarseness of the grid that can be tolerated. The following discussion 
will consider some of the aspects of the effect of grid size on the computer 
output, after noting some pertinent functional aspects revealed in the course 
of solving the problem. 

Figure 19 presents the geometry of a test case used to determine the 
effect of grid size on the results of the numerical simulations. The simu¬ 
lation was performed on UNIVAC to determine if a significant cost differential 
existed between it and the CDC series of machines. None was established, but 
differences in peripheral equipment were found to give each an advantage in 
different situations. This observation does underscore one feature that should 
be noted with regard to the use of machines from several different manufac¬ 
turers. The best way to handle the computer output is on magnetic tape. The 
problem arises with the variation in tape format and number of bytes per word. 

If differences do occur, conversion of the tape to the new machine format is 
required. This procedure is costly and, hence should be avoided by careful 
planning of computer usage. 

The main purpose of the simulation (Fig. 19) was to compare the effect 
of cell size variation on the basis of the situation presented in Fig. 18. 

A 5-cell-wide entrance was used in the simulation presented in Fig. 18 and 
a 10-cell-wide entrance in the other simulation presented in Fig. 19. The 
axial and radial velocities and the pressures at 10-msec time intervals at 
the locations noted in Fig. 18 are presented for both simulations in Table 6. 
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Fig. 19. Grid for Comparison of Variable-Cell-Size Procedure 
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Table 6 

COMPARISON OF THE EFFECT OF GRID COARSENESS 
ON THE OUTPUT OF NUMERICAL SIMULATIONS 




Axial Velocity 

(fps) 

Radial Velocity (fps) 

Overpressure 

(Psi) 

--Pt 

Time 

(msec) 

5 cell 

10 cell 

5c/10c 

5 cell 

10 cell 

m 

5 cell 

10 cell 

5c/10c 

■ 

n 

464 

463 

1.002 

-11.5 

- 7.4 

1.554 

4.95 

4.98 

.993 

WM 

MM 

610 

642 

.950 

13.3 

- 21.4 

- .621 

3.11 

3.31 

.940 

■ 

wm 

608 

616 

.987 

17.8 

- 22.6 

- .7876 

3.91 

4.41 

.887 


40 

432 

394 

1.0964 

13.2 

- 18.2 

- .7253 

7.45 

8.79 

.848 


10 

484 

515 

.940 

-46.3 

-47.4 

.977 

3.42 

3.42 

1.000 

o 

20 

541 

590 

.917 

-83.5 

-118. 

.708 

1.88 

2.00 

.803 

A 

30 

503 

509 

.988 

-83.9 

-112.7 

.744 

2.46 

3.05 

.807 


40 

239 

184 

1.2989 

-44.7 

- 53.8 

.831 

7.03 

8.27 

.850 


10 

243 

235 

1.034 

19.4 

11.8 

1.644 

2.53 

2.71 

.934 

o 

20 

514 

491 

1.047 

22.5 

17.6 

1.278 

2.61 

3.97 

.657 

«3 

30 

574 

601 

.955 

16.0 

17.2 

.930 

4.03 

3.98 

1.013 


40 

602 

676 

.891 

3.4 

6.01 

.566 

6.82 

6.48 

1.052 


10 

232 

219 

1.059 

73.4 

72.5 

1.012 

2.54 

2.63 

.966 

A 

20 

448 

441 

1.016 

87.3 

105 

.831 

2.61 

2*79 

.935 

*± 

30 

449 

482 

.932 

66.8 

112 

.596 

4.00 

3.79 

1.055 


40 

418 

439 

.952 

4.08 

12.8 

.3188 

6.54 

6.12 

1.069 


10 

62.6 

58.9 

1.063 

3.53 

1.5 

- 

1.05 

.93 

1. 129 

c 

20 

176 

177 

.994 

3.68 

2.38 

- 

3.64 

3.79 

.960 

0 

30 

128 

123 

1.041 

12 . 

10.4 

1.154 

7.56 

7.88 

.959 


40 

292 

284 

1.028 

11.5 

13.5 

.852 

7.88 

7.92 

.995 


10 

52.9 

51.9 

1.019 

10.1 

9.6 

1.256 

.93 

.84 

1. 107 

a 

20 

170 

171 

.994 

19.6 

17.1 

1.146 

3.77 

3.87 

.974 

t> 

30 

102 

105 

.971 

74.1 

70.8 

1.0466 

7.42 

7.65 

.970 


40 

223 

230 

.970 

71.6 

86.1 

.832 

7.83 

7.83 

1.000 


10 

10 

9.47 


- .0065 

1.000 

0 

0 

1.000 

n 

20 

10.4 

11.7 


4.85 

4.06 

1 .,195 

6,96 

4. 19 

.945 

/ 

30 

22.3 

18.1 


7. 

7.03 

.996 

7.60 

7.82 

.972 


40 

45.1 

39. 


10.6 

10 . 

1.06 

8.85 

8.67 

1.021 


10 

- .71 

.755 


- 1.53 



0 

0 

1.000 

8 

20 

10.3 

11.8 


30.6 

26.3 

1.163 

4.02 

4.11 

.978 

30 

19.0 

15.2 


53.4 

53.3 

1.002 

7.58 

7.77 

.976 


40 

32.3 

33.0 


66 .5 

66.8 

.996 

8.34 

8.28 

1.007 
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A review of the data presented in Table 6 indicates some regions in which 
deviations do exist between the 5- and 10-cell-wide entrances. In the entrance 
region (points 1, 2, 3) the inflow in the 5-cell-wide entrance persists for a 
slightly longer period than for the 10-cell case. One explanation for the 
longer flow is the poor representation of rarefaction waves by the coarser 
grids. 

A more comprehensive method of reviewing the data would be to consider 
the mass and energy influxes into the chamber. These basic parameters deter¬ 
mine the chamber conditions; hence by comparing the influx across the entry 
at comparable times, a comparison of the effects of the variation in coarse¬ 
ness of the mesh in the entrance can be made. The values of these parameters 
can be determined from the thermodynamic and flow parameters in the code out¬ 
put by summing the values determined for each cell. 

The values of mass and energy influx determined for any time step are, 
in effect, instantaneous values. Accordingly, the computed influxes are mass 
and energy rates at that instant, as opposed to quantitative values. Quanti¬ 
tative values were determined by combining the mean influx rate over the inter¬ 
val with the influx area and the time increment between intervals. This gross 
averaging process was sufficiently accurate, because the influx curves could 
be accurately approximated by linear segments for each computational interval. 

The values of mass rate were computed using the product of the vector 
velocity and the density. The energy influx was computed using the energy 
equation in the form of kinetic and potential energies (as discussed in Ref. 2) . 
The results of these computations are plotted in Figs. 20, 21, 22, and 23. The 
cumulative influxes defined in these curves are then the total influx of mass 
or energy to that instant in time. 

The simulation with the 10-cell entry produced higher influxes of mass 
and energy. The influxes of mass and total energy and potential energy had 
a maximum variation of approximately 3,5 percent between the two simulations, 
whereas the variation in kinetic energy influx was approximately 20 percent. 
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Fig. 20. Cumulative Mass Influx as a Function of Time 
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Fig. 21. Cumulative Total Energy Influx as a Function of Time 
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Fig. 23. Cumulative Kinetic Energy Influx as a Function of Time 
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The percentages would indicate that the variation in predicted pressure pro¬ 
files in the simulated room were not drastically affected by the number of 
entry cells. The kinetic energy variation was more pronounced, and velocity 
distributions and profiles would accordingly be expected to. vary. 

Figures 20 to 23 provide cumulative influxes as a function of time. 
Variations in the influxes of energy and mass can occur in large steps in 
short time intervals making it difficult to establish the error source. The 
dilemma is caused by the errors being carried forward cumulatively. To pro¬ 
vide a better insight into the process, the mass and energy influxes in 
specific computation time intervals are presented in Figs. 24 through 27. 

The largest deviations in the plotted data occurred in the kinetic energy 

influxes in the 15 to 20 msec time interval. The 10-cell simulation had a 
2 

37.7 x 10 ft-lb kinetic energy influx and a 0,0189 lbm mass influx, as com- 
2 

pared to 26.7 x 10 ft-lb and 0.0175 lbm influxes for the 5-cell simulation. 
These influxes correspond to average inflow velocities of 632 fps and 552 fps, 
respectively, with corresponding dynamic pressures of 3.32 psi and 2.68 psi. 

The variation in dynamic pressure corresponds to about 20 percent at this 
peak interval. The variation for most of the time range would be considerably 
lower. 

One interesting feature is that the kinetic energy, although a small 
portion of the total energy involved, does play a large role in determining 
the drag forces of the flow. A brief review of why the discrepancy occurred 
is accordingly appropriate. The basic reason for the variation is the in¬ 
ability of the coarser grid to accurately reproduce gradients in the flow. 

These gradients of interest are the entering shock wave and its subse¬ 
quent reflections, the rarefaction waves propagating around the entry and 

through the chamber, and the rapidly curving streamlines of the inflow. The 
first two factors cause delays in the changes in the flow variables, but 
should not form long-term variations. The problems cannot, however, be com¬ 
pletely isolated because the inflow is dependent upon the movement of these 
waves. The waves initiate the strong inflow that develops the strong gradients 
in the streamline curvature. 
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Fig. 25. Total Energy Influx as a Function of Computation Intervals 


3-58 








KINETIC ENERGY INFLUX IN EACH TIME INTERVAL (10 2 ft-lb) 





UR5B 755-3 


The problem of the rapid change in curvature can be considered by review¬ 
ing the contraction equations for a vena contracta. To arrive at an estimate 
of the curvature of a streamline, the maximum curvature for a free jet of 
incompressible fluid is selected. The maximum curvature would occur on the 
free streamline. Conformal mapping via the Schwarz-Christofel theorem as 
applied by Chaplygin will provide a functional form for the streamline. Dis¬ 
cussions of these methods can be found in Refs. 45 to 53. 


The parametric equations of the steady free surface, first developed by 
Lamb, are given by Ref, 45 as 


x 

L 


7r 

2 + 7r 


2 

2 + 7T 


cos 6 


(14) 


y 

L 


2 

2 + 7T 



sec 6 + tan 0 ) - sin 0 


(15 


\ 

/ 


where the terms in the equations are defined in Fig. 28. Define a cell width 
w equal to L/n^, where n^ equals the number of cells in distance L, and assume 
the cells are square. Since these equations cannot be worked with easily, a 
graphic approach to the interacting effects of vena contracta and cell size 
will be used. 

The variation in the cell size as shown in Fig. 28 affects the coverage 
in the first cells. For the three cases shown, w = O.lx/L, 0.2x/L, and 
0.25x/L, i.e., corresponding to 10, 5 and 4 cells per half-width of entry, 
none of the cases provided a capability to adequately represent the flow. The 
value of the flow calculated for each cell is at the center of the cell. In 
the cases shown in Fig. 28, the centers of all the cells are outside of the 
vena contracta, in a zone in which there is no flow. 

The averaging process by which the conditions of the cell boundaries are 
calculated (see discussion in Section II and Appendix A) can provide a method 
for accounting for this type of variation. Due to the rapid change in flow 
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properties and the rapid curvature of the streamline, the flow field cannot 
be accurately represented by the cells shown. The reason for the greater flows 
in the 10-cell case can be deduced from Fig. 28. 

In the smaller cell size (n = 10) , the first cell against the barrier 

c 

does not contain the free streamline; but a cell below it would contain 85 
percent of the corresponding flow field for that zone. The first cell area 
contained 25 percent of the flow inside the vena contracta in that cell. For 
the case of = 5, approximately 30 percent of the first cell's area is in the 
flow zone of the vena contracta. This small percentage corresponds to an 
inability of the larger single cell to represent the flow accurately. In the 
comparisons of energy and mass influx, the result was a higher influx of mass 
and energy for the 10-cell case, which could more accurately represent the 
inflow. The 5-cell case represented a greater spreading of the inflow to lower 
values. 

The reduced influx in the 5-cell case appears to correlate more closely 
to experiment. The reduced influx corresponds to some of the energy-dissipation 
mechanisms internal to the flow fields which the numerical simulations do not 
account for. The dissipation mechanisms appear to be lessened in larger scale 
shelters, but sufficient experimental evidence is not available to completely 
document this observation. 

In that the deviation between 5- and 10-cell data seems to favor use of 
the smaller number of zones, it would seem possible for a coarse-grid code 
to be used to predict flows in shelters. Additional analytical and experi¬ 
mental work should be undertaken to fully ascertain the implications. If the 
present indications are supported, a fairly simple and very general two-dimen¬ 
sional code could be developed as a general prediction tool for the shelter 
designer. 
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Section 4 
FLOW PROCESSES 


This section will review several aspects of the general flow process 
prior to discussing a qualitative model of the flow. Based on prior work 
(Ref. 2), the following general model of the total inflow process is pre¬ 
sented as preliminary groundwork for the individual flow processes. The 
qualitative or conceptual model will then be formulated using these results, 
the numerical simulations, and implications drawn from experimental studies. 

The airblast wave, on reaching the exterior of the shelter chamber, es¬ 
tablishes a high pressure reservoir which provides the energy for driving the 
inflow of gas into the shelter. The airblast initially enters as a shock wave 
which propagates and expands through the chamber. This transient process is 
seen as a perturbation on the long term filling process of the chamber. 

The inflow is formed by an adiabatic expansion from the reservoir to the 
chamber conditions. The influent gas is a high-velocity, highly directed flow 
which is in effect a transient jet. The inflow is the driving force for a 
long term circulatory flow field. The pattern is altered if the chamber is 
long relative to the side in which the entrance exists. For this case the flow 
approaches that found in a tunnel in which the circulatory flows are replaced 
by axial flows. 

THE TRANSIENT INFLOW PROCESS 

The initial expansion of gas into the room is a two-part phenomenon. The 
first is the transmission of the incident shock wave into the room and its 
subsequent expansion. The second part is the evolvement of the flow field 
behind the shock wave. In early times (in the range of t = 0 to t = d/a 2 
where d is the slit width and a is the sound velocity behind the incident 
shock wave) the transmission process governs, but with time decreases in 
effect until the flow expansion behind the shock wave governs the flow. In 
this section some of the characteristics of the formation process will be 
reviewed. 
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A specific precise mathematical solution for the unsteady flow of a jet 
of compressible fluid could not be found in the literature. The closest 
approximation to the process was found in a presentation by Curie (Refs. 52 
and 53) on the expansion of an inviscid incompressible fluid from a vessel. 
Although an error occurs for the compressible flows being considered, this 
work will be reviewed and used to provide further inputs to the conceptual 
model being developed. 

The unsteady flow representation is based on the classical Helmholtz 
conformal mapping procedures found in Refs. 45 through 51. These works are 
the outgrowth of Lamb’s now classical solution (Ref. 51) to the problem of 
defining the free flow of an inviscid and incompressible flow through a slit 
from a large reservoir into free space. Although the flows are incompressible, 
the close correlation of the incompressible flows to compressible flows can be 
found in the work based on Chaplygin’s techniques (Ref. 50). The steady state 
representation, such as the one shown in the preceding section, is nnt directly 
applicable to the transient portion of the compressible flows encountered in 
shelter flows. 

The unsteady flow analysis developed by Curie (Refs. 52 and 53) does 
provide some interesting, though limited, insight into the transient flow 
definition problem. A brief review of his work will accordingly be undertaken 
and its significance relative to the conceptual model will be discussed. 

The basic analytical procedure in conformal mapping is the transformation 
of the physical plane to the hodograph plane using a complex velocity. The 
rationale for this procedure is that 9 (the flow vector direction) is then a 
constant along all solid boundaries and q (the magnitude of the flow vector) 
is constant on free streamlines. The procedure, a further description of the 
mapping process, and examples of applications of the theory can be found in 
the references. 
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Curie's contribution extends the steady state analysis of jets to include 
their transient formation. The solutions apply to the expansion of an incom¬ 
pressible fluid into free space. Since the purpose of this study is to gain 
insight into the processes, the use of the incompressible approach seems jus¬ 
tified because analystical work on compressible flows is not available. It must 
be stressed that the approach will only provide analytical approximations, very 
much qualitative as opposed to quantitative. 

, -At 

Curie's method is to expand the velocity potential 0 in powers of e 
in the form 

0 (x,y,t) = 0 o (x,y) + e - * 1 0 1 (x,y) + 0 (e (16) 

Equation (16) is basically the pertubation of the steady state solution 
(x,y) by time-dependent terms including an Eigenvalue X which is determined 
by the boundary condition at infinity. Curie also defined the spatial vari¬ 
ation of the unsteady flow [£> (s,t)] from the steady state boundary by the 
expression 

0 (s,t) = 6 (s) e ^ + 0 (e 2 ^) ( 17 ) 

Through these expressions Curie has altered the physical boundaries of the 
flow to correspond to the transient situation. The mathematics involved in 
the process can be found in Refs. 52 and 53. 

The first result Curie predicts in his work is the existence of an 
M-shaped discharge in the jet, shown in the upper portion of Fig. 29. The 
strong influx along the edges of the entry provides the reverse effects of 
the vena contracta which eventually forms with its strong flow in the central 
region. Curie points out that extension of this analysis would require a 
more detailed knowledge of the upstream conditions. 
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Fig. 29. Curie's Representations of Early Inflow Process 
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The lower portion of Fig. 29 shows the outline of the jet shape at vary¬ 
ing times after inception of the flow into the downstream region. The dashed 
curves show improved estimates of the curvature of the nose of the jet which 
Curie points out is inversely proportional to the square of its velocity and 
that the limiting radius of curvature approaches four tenths of one-half the 
slit width. The lower dashed curve in the bottom portion of Fig. 29 represents 
the shape of the jet when the nose has propagated more than a slit width down¬ 
stream. 

The conditions of the transient compressible flow situation negate many 
of the results reached in the preceding analysis. The initiation of the flow 
by the propagation of a shock provides a downstream flow field for the jet 
to flow into. Similarly, the reflection processes in the upstream region 
also contribute to the formation of the inflow process. Unlike the incompres¬ 
sible solution in which free space was assumed, the downstream region contains 
mass which must be moved by the influent gas. 

The above, although negating the possibility of a direct application of 
Curie’s work to the shelter inflow situation, does indicate features in the 
flow which may be of interest. These features are the result of the flow 
dynamics and should occur in some form in the compressible flow. Of particular 
interest is the profile of the influx as it moves into the chamber. The 
narrowing of the inflow in these early stages corresponds to the lack of spread 
of the jet observed in the large scale experimental studies discussed in 
Appendix C. The profiles of the flows in the numerical simulations resemble 
the described phenomena, but superpositioned on the shock induced flow field. 

Further work in this area was not undertaken. Future studies should re¬ 
view this area, for if the effect predicted by Curie does exist in shelter 
flows, a limitation in the spreading of the inflow would result during the 
transient flows of the filling process. The following review of the transition 
processes considers the compressible gas effects via a one-dimensional approx¬ 
imation. 
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The inflow process in its initial stages for a compressible gas is 
comprised of a complex system of interacting waves. The process was de¬ 
scribed in Ref. 2 for a simple two dimensional geometry using the wave 
systems which form upon the incidence of the shock wave at the tunnel en¬ 
trance. Figure 30 presents the wave configurations for both a narrow and a 
long entry. To gain further insight into the transient inflow process, the 
two-dimensional system will be transformed into a one-dimensional system which 
approximates the described wave systems. 

For the purpose of this discussion, the entry will be extended through 
the chamber in the form of a long tunnel. The extension is in effect a model 
of the jet that forms without any components interacting with the remainder 
of the chamber. The reservoir conditions will also be developed on a one¬ 
dimensional basis with constant pressure conditions imposed to form the source. 
The resultant geometric and wave diagram schematics are presented in Fig. 31. 

The increase in the velocity of the reflected shock wave at early times 
corresponds to the closure of the reflected wave and the gradual approach to 
equilibrium conditions. The wave diagram shown in Fig. 31 is for an infinite 
reservoir. For a finite reservoir, the rarefaction fan shown would have a 
continuing interaction with the reflected wave until equilibrium conditions 
were reached. The equilibrium conditions would result when the upstream mass 
influx into the system was balanced with the mass influx into the chamber 
through the entry. 

A decrease in the velocity of the transmitted shock wave begins early 
in the transmission process. The decrease corresponds to the interaction of 
the rarefaction wave with the transmitted shock wave (labeled 1 and 5 in 
Fig. 30, respectively). The rarefaction fan shown in the wave diagram of 
Fig, 31 represents the acceleration of the fluid from the reservoir into the 
jet corresponding to the pressure differential between the two regions. 
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(1) Transmitted incident shock wave 

(2) Reflected incident shock wave 

(3) Rarefaction wave from exterior corner 

(4) Reflected wave of Mach configuration 

(5) Rarefaction wave from interior corner 

(6) Regular reflection wave from wall 
interaction 

(Abstracted from Ref. 3) 




(Abstracted from Ref. 2) 


Fig. 30. Blast Transmission into a Narrow Entrance and the 
Equivalent Process in a Long Entry 
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31. A Simple Diagram of the Initiation of a Quasi-One- 
Dimens ional Process 
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The expansion process during the transitory period creates an interesting 
anomaly. If a gas is expanded from the reservoir conditons to the pressure 
in the jet (or behind the shock front), a difference in velocities results. 
Since two masses cannot occupy the same space, the influent mass acts as a 
piston on the gas in the chamber (or brought in by the incident transmitted 
shock wave). The contact surface is then the boundary of the piston from 
which emanates a series of compression waves. As the piston is a gas and 
not a solid, a series of backward-facing rarefaction waves must also emanate 
from the contact surface. 

In a steady flow these waves are absent as the shocked gas ahead of the 
contact surface has moved away and a steady expansion and flow condition has 
been reached. Under normal two-dimensional steady flow conditions, the fluid 
in the jet is a combination of influent gas and chamber gas, with the condition 
that at any cross-section the momentum is approximately constant. In the tran¬ 
sient chamber-filling situation this condition is never precisely reached, and 
only at certain times in the filling process is it even approached. 

The mathematical description of the inflow process is complicated by its 
nonlinearity. A quasi-steady approximation of the later portions of the 
inflow process was presented in Ref. 3 to provide some bounds on the process. 

The transiency of the flow, particularly at early times in the inflow, produces 
a wide variety of possible solutions depending upon the intermediate conditions. 
Several of these intermediate situations were approximated by quasi-steady 
flow analogies and were investigated in Ref. 2 and shown to provide valid 
local solutions. 

To illustrate the process on a transient basis, the following example 
was selected. The incident wave impinges normally upon an infinte reflecting 
wall containing the chamber entrance. The wave will be assumed to be fully 
reflected in the reservoir region and the incident wave will be transmitted 
into the chamber without any attenuation (the most unfavorable condition for 
this argument). A pressure differential across the entry exists corresponding 
to the difference between incident and reflected conditions. 
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Rather than develop similar derivations as found in Ref. 2, the values 
presented therein will be used and are presented in Table 7 for the preceding 
conditions. For the most favorable, or lower bound, condition, the variation 
between inflow and shocked flow values range from multiples of 2 to 4. The 
upper bound flows, a full adiabatic expansion, have multiples ranging from 
approximately 3 to 7, The piston effect must exist if mass is to be conserved. 


Table 7 

COMPARISON OF FLOW PARTICLE VELOCITIES 


Overpressure (psi) 
Incident Reflected 

u 2 (fps) 

* 

u 

e 

** 

u 

e 

1 

2.06 

53.4 

362 

195 

2 

4.23 

104 

502 

333 

3 

6.51 

152 

610 

416 

4 

8.90 

198 

730 

468 

5 

11.4 

242 

831 

503 

6 

14.0 

284 

921 

528 


u^ = flow velocity behind incident wave, assumed to be the 
same as transmitted wave 

u = expanded velocity for a pressure differential between 
incident and reflected conditions 

* = pre-contraction, refers to flow value computed at the 
outer entrance of the chamber (Ref. 2) 

** == max-contraction, refers to the velocity obtained in a 
free adiabatic expansion (Refs, 2 and 3) 

To illustrate the piston effect of the flow, we first consider a simli- 
fied case in the quasi-one-dimensional system. The inflow process must be 
modeled in the transformed one-dimensional coordinates for the more general 
case, hence the simplified case will use a start-up model, shown in Fig. 32. 
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In the preliminary model (shown in Fig. 32) an artificial system is estab¬ 
lished. To account for the time span in originating the transmitted and re¬ 
flected waves one dimensionally, in comparison to their two dimensional 
counterparts, a 0.4 x/d spread was assumed (a close approximation to the closing 
model given in Ref. 2). As the flow begins, this standing wave remains and 
defines the expansion between the reservoir and the chamber. In the two-dimen¬ 
sional case, as the wave systems formed, the transmitted shock wave (T) was 
weakened by its expansion in the chamber and the hole in the reflected shock 
(S) closed. The compression wave (C) and the rarefaction wave (R ) are the 
conditions imposed to account for the beginning of the process. These would 
be impressed in steps, over a period t = 0 to a i t/d=l, such as the one shown 
at t 3 . 

If we review Fig. 32, during a period prior to t ( such as t ) the follow- 

» 1 

ing sequence is observed: 

• An inflow corresponding to a 5,0 psi overpressure shock wave flows 
through the entrance and the reflected shock (S) is raised to an 
overpressure of 5.5 psi. 

• The inflow is then expanded through the standing rarefaction wave 
to the 5 psi overpressure behind the transmitted wave. 

Not shown in the wave diagram is the effect of the expansion on the flow behind 
the transmitted shock, i.e., the location of the contact surface. More specifi¬ 
cally, if the expansion is adiabatic, does the expanded flow velocity equal 
the flow velocity behind the transmitted shock wave? This condition must be 
met if a contact surface is to exist (across which flow velocity and pressure 
are equal) . 

Using the above properties of the flow, the conditions behind the reflected 
shock can easily be determined by transforming the coordinate system so as to 
stop the wave in laboratory coordinates and then writing the continuity equation 
across the wave. For the weak waves being considered, adiabatic conditions 
can be assumed, and using a to represent sound speed, p pressure, u particle 
velocity, p density, and y the adiabatic gas constant, the following equations 
can be written. 
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(18) 


(19) 


The expansion from condition 7 to 2 is adiabatic and hence not path-dependent. 
The energy equation for steady flow can be written as 


_y 

ri 


> • ?*(*) 


( 20 ) 


and as p 0 = p_ then for a completely adiabatic process a = a . Using the 

8 2* o A 

above relationships and combining them with the preceding equations, we can 
write o 


u - a 
8 8 


H#] 


- 5 (a 2 “ a 2 ) 

7 8 7 


( 21 ) 


Substituting values into Eq. (21) for the 5 psi incident overpressure condition, 
an expanded velocity of 313.5 fps is determined [the shocked flow velocity in 
region (2) is 241.5 fps]. The expanded flow is going faster than the shocked 
flow, hence it is acting like a piston. 


To follow this line of reasoning a wave diagram will be constructed (Fig. 
33) showing the initial sectors of the inflow formation process. The quasi- 
one-dimensional construct will again be used to provide a model for the process. 
The rapidly increasing complexity of the wave diagram precludes a full solu¬ 
tion for the full period of inflow into the chamber. The preliminary phases 
of the wave diagram do, however, provide valuable insight into formation of 
the flow process. 


A slight modification to the standing rarefaction wave has been made. 

As shown in the wave diagram, the standing rarefaction waves follow particle 
paths (zones 2', 3, 4, 5, and 6). An artificial separation of the variables 
occurs in zones 2 f , 2 and 12 and again between 13 and 16. A series of 
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Fig. 33. Wave Diagram Approximation for Quasi-One-Dimensional Inflow Simulation 
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compression waves would emanate from gas particles expanding into zones 2 and 
12 from 2 1 or 16 from 13. Rather than use a series of waves a single wave was 
used to represent these waves, hence an inconsistency occurs at the boundaries 
mentioned. 

Tables 8 shows the values of the flow, state, and Reinman invariants for 
each zone. Adiabatic conditions were assumed across all waves, introducing a 
small error but simplifying the calculation procedure. In light of the assump 
tions made to form the quasi-one-dimensional process, the small error attribut 
able to entropy increases is negligible. This is in keeping with the purpose 
of the exercise, that of providing insight into the flow process. 

The flow (u) across the expansion (regions 7, 6, 5, 4, 3, 2, 2* ) is seen 
to increase above the value for region 2 as was indicated by Eqs. (18) to (21) 
The result is a ^^compression wave moving into region 12. The pressure and 
velocity in regions 13 and 14 were matched by adjusting the value of the Q 
wave moving into regions 2 f and 3. The expansion from 4 to 3 creates a 
velocity (u) whose corresponding pressure is greater than the pressure in 
either region 3 or 13 and, therefore ^and Q compression waves form. 

The expansion of the shock wave in the chamber is simulated by the 
rarefaction waves between regions 2 and 12 and 14 and 16. The closure of the 
reflected shock wave by compression waves occurs between regions 7 and 8 and 9 
and 10. In the investigation, the process discussed was varied by assuming 
different conditions at various interfaces. The results obtained from the 
other conceptualizations of the flow forming process did not vary appreciably. 

Arguments can be produced to fault the assumptions of the preceding ana¬ 
lysis. The results presented were for what appeared to be the most realistic 
of the wave diagrams constructed. Some similar diagrams for a duct exiting 
into free air can be found in Ref. 54. These wave diagrams carry out the 
computations to later times and show the pulsating nature of this form of flow 
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Point 

r 

(5a+u) 

a 

(5a-u) 

u 

a 

u+a 

u-a 

p 

(psi) 

1 

5 

5 

0 

1.0000 

1.0000 

- 1.0000 

0 

2 

5.42941 

5.00179 

0.213810 

1.043120 

1.25693 

-0,82931 

5.0 

3 

5.4829808 

4.9563192 

0.2633308 

1.04393 

1.3072608 

-0.7805992 

5.1 

4 

5.47140 

4.975439 

0.2479806 

1.044684 

1.29266 

-0.796703 

5.2 

5 

5.458826 

4.995474 

0.231676 

1.04543 

1.277106 

-0.81375 

5.3 

6 

5.44524 

5.0164 

0.21439 

1.04617 

1.26056 

-0.83178 

5.4 

7 

5.42983 

5.03937 

0.19523 

1.04692 

1.24215 

-0.85169 

5.5 

8 

5.42941 

5.00179 

0.213810 

1.043120 

1.25693 

-0.82931 

5.0 

9 

5.42983 

5.075986 

0.1769135 

1.05058 

1.227494 

-0.873667 

6.0 

10 

5.42941 

5.061828 

0.18379058 

1.0491239 

1.2329144 

-0.8653333 

5.8 

11 

5.410396 

5.00179 

0.2043028 

1.0412186 

1.2455214 

-0.836916 

4.75 

12 

5.391170 

5.00179 

0.194690 

1.039296 

1.233987 

-0.844606 

4.5 

13 

5.4829808 

5.00179 

0.240595 

1.04847708 

1.28907208 

-0.807882 

5.7 

14 

5.4829808 

5.00179 

0.240595 

1.04847708 

1.28907208 

-0.807882 

5.7 

15 

5.4638608 

5.00179 

0.2310354 

1.0465651 

1.27760048 

-0.81552968 

5.18 

16 

5.4447408 

5.00179 

0.2214754 

1.04465308 

1.26612848 

-0.82317768 

4.93 

17 

- 







18 

- 







19 

5.47140 

5.00179 

0.234805 

1.047319 

1.282124 

-0.8125140 

5.54 

20 

5.47140 

5.00179 

0.234805 

1.047319 

1.282124 

-0.8125140 

5.54 
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velocity. The process changes the shape of the jet's velocity profile through¬ 
out a transition region which ends with the evolvement of velocity profiles 
that will be maintained for the remainder of the flow. The balance of the 
acceleration of the outer gas and the deceleration of the inner gas can be 
considered to be a redistribution of momentum through a mixing process. 

Variations of these properties of basic jet flow can be established, but 
they were not considered in Ref. 3. The jets were further categorized as in¬ 
compressible, turbulent, free, viscous, subsonic, and submerged. In the 
present discussion only the plane two-dimensional jet will be considered. 
Well-documented and relatively simple solutions exist for incompressible jet 
flow with which also have all of the other characteristics listed above; 
hence the effect of incompressibility is of particular importance. 

The magnitude of the effects of compressibility and heat transfer on the 
flow of subsonic jets were found to have been investigated by Abramovich 
(Ref. 55) . The result of his investigations indicated that if heat transfer 
was low; i.e., the temperature difference between the core and the stationary 
gas was small, the effect on the flow's velocity profile was small. Similarly, 
he determined that the effects of compressibility were yet more limited for 
flows below Mach 1. The work in Ref. 2, showed these conditions are approxi¬ 
mated by shelter flows. 

The centerline velocity distribution for a plane jet is presented for 
several driving pressure differentials in Fig. 35. These centerline distribu¬ 
tions were coupled to a full representation of the flow field in Fig. 36 for 
a pressure differential of 12 psi across the entranceway. A striking feature 
of these flows is the distance the jet flows before there is an appreciable 
decrease in the velocity; for example, the core extends into the chamber 5.5 
diameters. The transverse growth of the jet is more restricted and the intense 
flows extend only 1.5 to 2.0 diameters from the centerline of the chamber. 


4-20 



AXIAL VELOCITY (FPS) 


URE ■ 


755-3 



Fig. 35. Centerline Velocity Distribution 
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Reference 3 also discussed some general characteristics of jets. Since 
the pressures in the free jet are relatively constant (Ref, 56) and the tem¬ 
perature variations small, the density of the jet flow does not vary greatly. 
Transverse flows occur only as perturbations in the flow field and present no 
lasting effect. A time period elapses between entrance of the shock wave into 
the chamber and the establishment of the fully developed jet as was discussed 
earlier in this report. Experimental evidence (Refs. 3 and 4) showed the 
penetration of the fully formed core into the chamber ranges from 3.5 to 5.5 
effective entrance widths. The term "effective width refers to the contraction 
of the inflow in the entrance corresponding to the formation of a vena contracta. 
The experimentally measured contraction ratios appear to correspond to values 
established by steady flows. 

The preceding synopsis was a condensation of the initial estimate of the 
inflow process into a two-dimensional shelter space in the first of the present 
OCD study series (Ref. 3). Recent work has shown variations in the inflow 
process from the initial model of the flow do exist and warrant attention 
before considering an improved model of the inflow process. Earlier in this 
section, analytical models of the inflow process were discussed that showed 
the dynamics of the inflow process to play a stronger role than was assumed 
in the first work. The interaction of the gas ahead of the accelerating 
influent gas and the continuing interaction between the components of the flow 
field were shown to play a major role in the evolvement of a steady flow jet. 

Another variation from the quasi-steady characterization is that a flow 
field does exist within the chamber which reduces some of the transverse 
velocity discontinuities present in the steady flow free jet. The existence 
of the internal flow fields and their magnitudes are represented in the data 
discussed in Section 3 and in Appendix C. These flow fields are related to 
the piston effect (discussed earlier in this section) in which the motion of 
the gas ahead of the expanded inflow is driven by it; forming as one boundary 
condition on the inflow, the expanding shock wave in the chamber. The propa¬ 
gating shock wave will not terminate at the steady-state jet f s shear flow 
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boundary, but rather will spread across the chamber with decreasing strength 
as documented in Refs. 3 and 4. Accordingly, the quiescent gas surrounding 
the inflow in the shelter situation does not exist in the pure form of a steady 
state free jet. 

The preceding comments stress that the dynamics of the early portion of 
the inflow process play a significant role in forming the flow field. Curie T s 
work provides an approximation of the form of the early inflow process. This 
description indicates the dynamics of the transient inflow may inhibit the 
spreading of the influent jet of gas. These analytic indications, when coupled 
with the results from the large scale experiments indicate that a further 
review of the spreading process should be undertaken. These conclusions are, 
however, contraindicated by the data from the scale model experiments which 
show a pronounced turbulent spreading of the inflow. The opposing results 
stress that a better understanding of the shelter flow process is needed. 

To better understand the shelter inflow, a literature review was under¬ 
taken. Information was sought in several general areas: 

• Characteristics of flows with the general character of 
shelter flows 

• Shear flow 

• Turbulence, mixing, and transition processes 

• Impinging jets and wall jets 

The literature review was then edited and the results presented in the 
bibliography in Section 7. 

A review of the bibliography provides an insight into the complexity 
of the flow processes involved. Not only the volume of work, but also that 
it primarily represents quasi-steady or steady flows establishes the diffi¬ 
culty in gaining an understanding of the flow sufficient to develop a predic¬ 
tive model. Based on the literature review, several different approaches to 
studying the problems associated with the shelter flows were undertaken. 


i 
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Due to the limitations and time contraints on the project, these investiga¬ 
tions cannot be reported, but rather the combined overview gained will be 
presented and one very general approach to several major problem areas will 
be defined. 

The basic hypothesis for the shelter inflow model is that a jet of gas 
flows into the shelter behind the expansion of the incident shock wave through 
the chamber. In the present argument, it is assumed that the surrounding 
air is stationary and that the jet is a potential core of laminar constant 
velocity flow. This situation creates a maximum velocity discontinuity which 
will be considered to be an infinitely thin line of separation, a laminar 
shear layer. The situation is then in the form of the classical "Helmholtz 
instability" problem. 

For the divided flow, the shear layer can be considered to be an infinitely 
thin vortex sheet. The vortex sheet represents the introduction of an insta¬ 
bility which caused the shear layer to grow to finite proportions. The rate 
of growth of the instability corresponds to a decrease in velocity in the 
potential flow and an increase in the velocity of the stationary gas. The 
conditions just described are similar to those described for the profile of 
a steady flow jet described earlier. 

Of particular interest to the shelter flows is the growth rate of the 
instability, i.e., how rapidly does the momentum transfer from influent to 
stationary gas occur. The slower the spread the less the deceleration of the 
influent gas and the less the spread of the flow profile normal to the inflow 
axis. Several approaches to this problem were undertaken, the simplest and 
most general of these was to consider the shear line as a surface from which 
a boundary layer evolved in both directions. In this format the evolvement 
of a boundary layer behind a shock wave on a flat plate can be used to 
approximate the flow processes. 

The flow at any point on the flat plate begins with passage of the shock 
wave over that point. The flow can be turned into a steady-state problem by 
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considering the point to be the leading edge of a plate inserted into the flow 
stream, as shown in Fig. 37. For the initial purposes of this discussion, 
the flow field behind the shock wave will be considered to be that induced 
by a step shock wave. The following analysis was carried out in Ref. 57 for 
scaling effects on supersonic flows over antennas and will be reapplied and 
expanded for the present problem. 

The displacement thickness of a boundary is given by the relationship 


<5*(x) 



( 22 ) 


z is a temperature and viscosity dependent constant 
V is the dynamic viscosity 
u is the speed of the fluid 

x is the distance from the leading edge of the plate 
6* is the displacement thickness of the boundary layer defined by the 
integral 



_ 

P U/ 

GO ' 


dy 


U is the free stream velocity 
p is the density 
p is the free stream density 

oc J 

y is the coordinate perpendicular to the plate 
0 is a constant 


Using the nomenclature in Fig. 37, the dimension (x) can be replaced 
by U t, and the flow stream properties (subscripted by 2) by the Rankine- 
Hugonoit conditions behind the incident shock wave, or 


<5*(t) = 



(23) 
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The above expression can be rewritten in terms of the thermodynamic and 
aerodynamic parameter ratios across the shock front as 


<5*(t) 



(24) 


where T 21 is the density ratio, U 21 the velocity ratio, and W X1 the shock 
front Mach number. The rate of displacement thickness growth is given by 


d 6*(t) 
dt 



(25) 


Both expressions show that the boundary layer profile does not scale linearly 
with time. 

The work in Ref. 57 extended the preceding expression to consider the 

nonscaling character of the boundary layer in terms of a model and prototype 

situation. A linear scale factor (f ) was defined by the ratio t /t . time 

s m p 

in the model (t ) and prototype (t ) reference frames, respectively. If the 

P 

initial conditions and the shock strengths are equal, then 


and 


<5*(t> nr 

m / m 

d’TtT = J T = /f s < 26 > 

p V p 


d <5*(t) /dt 

m 

d <5* (t) /dt 

P 



(27) 


To see what this would mean in terms of the boundary layer growth 
between a model and prototype situation, the BRL 4-in. shelter models (Fig. 2) 
were selected. The prototype room selected for comparison was 21.3 ft in 
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length. The linear geometric scale is 64 for this situation, which also is 
the factor for scaling time between model and prototype situations. Under 
these conditions the boundary layer in the modeled situation would grow 
eight times faster. 


The above result helps to explain why the pronounced turbulent spread 
observed in the model data was not observed in the full-scale test results. 
The indications are then that the width of the mixing zone is reduced, but 
as the momentum cross-section must remain constant, the flow must be more 
intense. These results appeared to be reproduced by the preliminary work 
undertaken using the spread of a Helmholtz instability as the basis for the 
boundary layer growth. Preliminary work using the momentum transfer argu¬ 
ments, based on either Prantl's mixing length theory or Taylor's vorticity 
transport theorem, appeared to be yielding supportive results when this study 
ended. 

The flows in the jet are quite intense, hence the boundary layer (which 
is being used to represent the jet's spread) will grow more quickly than for 
the flow conditions behind the shock wave. Since the disturbance is still 
initiated by a shock front, the shock front velocity will be retained, but 
the other properties will be taken from the conditions of an adiabatic expan¬ 
sion that forms the jet. Rather than use the complex expression described 
above, the following approximate formula can be applied (Ref. 58). 


1.75 



(28) 


and if x = U t , then 
s 


<5*(t) 


1.75 


i 


V V t 
s 


u 

OO 


(29) 
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If we use the assumptions that the jet is formed by an adiabatic 
expansion to the initial chamber pressure, and that the shock front velocity 
is approximately equal for all shock strengths (selecting 1300 fps which has 
a maximum error of approximately ±15 percent in the range of interest) then 
Eq. (29) can be written as 


<5*(t) ^ 0.804 Vt7u (30) 

QO 

The boundary layer thickness, 5(t), is approximately two and one-half times 
greater. Values for these parameters as a function of time and free-stream 
flow speed are presented in Fig. 38. 

» 

The curves in Fig. 38 support the indications of the nonscaling charac¬ 
ter of the flow. To ascertain the magnitude, the two sets of experimental 
data used earlier will be considered (the BRL 4 x 4 x 4-in. model with a 
1-in. opening and the URS shock tunnel test of a 12 x 15 x 8 ft room with a 
2.7-ft opening). The last comparison time used with the BRL model was 
approximately 750 jj,sec. On a linear time scaling basis, the comparable 
time in the large-scale model would have been 18 msec (using the shortest 
room dimension). In 750 (ji sec the boundary layer thickness would have ranged 
between 0.02 and 0.07 in. in thickness. If we assume that this growth approxi¬ 
mates the spread of the jet, then up to a 0.14 in. widening of the inflow for 
a 1-in. opening would occur. The full-scale room for an 18 msec fill time 
would show a maximum growth of 0.6 in. on a 27-in. door, i.e., the full-scale 
situation would exhibit a growth l/6th less with respect to the opening width. 

Linear scaling does not also apply to the filling of the chamber. The 
fill time (t^) in milliseconds can be approximated (Ref. 44) for a two- 
dimensional room by the expression 

V 

t = —2- 

f 2A 

e 


2W 


( 31 ) 
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Fig. 38. Approximated Boundary Layer Thickness Parameters for a Free-Flow Jet 
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where V = volume of the chamber of floor area A and A = area of the 
c f e 

entrance of width W (all dimensions are in feet). The fill time for the 

e 

4-in. model would then be 667 jjtsec and 38 msec for the full-size room (linear 
scaling on the basis of opening.size would give an 18 msec fill time). The 
argument of the nonscaling character of the flows in the models does not 
reduce their value in understanding the flow; it simply cautions against 
direct scaling and extrapolation to full-scale flows. 

The boundary layer arguments presented apply equally well to the flows 
along the shelter walls, ceiling, and floor. The inflow is therefore not a 
free jet, but rather is more nearly a wall jet. In Fig. 38 the indication is 
that the effect of the walls in the form of the boundary layer build-up would 
be minimal for the large scale, but would have observable effects in the model 
flows. In the review of the data undertaken earlier, a decrease in the 
accuracy of the smoke grid experiments was noted with time as was a widening 
of the smoke grids. These occurrences are felt to be partially due to the 
boundary layer growth. 

The growth of the boundary layer or the spread of the turbulence are 
energy dissipation and momentum transfer mechanisms in the flow. The jet 
flows into the shelter and to a pressure that is between the ambient and 
filled pressure values. The deceleration force, in terms of pressure gra¬ 
dients, is therefore not sufficient to reduce the flow in the chamber to 
zero. The momentum of the flow must therefore be transferred through mechan¬ 
isms such as turbulence. The reduced spread of turbulence in the large 
scale rooms results in a long-term circulatory flow well past the filling 
period of the room (as documented by the URS shock tunnel tests). 

There is a reduction in intensity as the inflow strikes an opposing wall 
and turns. The flow stream spreads and, since momentum must be conserved, a 
reduction in flow intensity occurs. The work unit ended as this portion of 
the flow was being reviewed. Some observations were noted from the prelimin¬ 
ary work and will be briefly outlined. 
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The turning of the flow at an opposing wall and the long-term circula¬ 
tion appeared to be a function not only of the dynamics of the flow deflec¬ 
tion but also of the wave dynamics within the chamber and the balance between 
kinetic and potential energy. The momentum generated by the influent air 
appeared to be retained with a large kinetic energy component well after the 
filling process was terminated. 

The moving pictures taken in the large-scale room tests in the URS 
shock tunnel corroborated the preceding results. A circulatory flow was 
observed at 500 msec, more than an order of magnitude in time past the filling 
period. This long circulatory period was partially induced by the cyclic 
reflections in the closed tunnel, but the reflection process was not solely 
responsible. For example, circulatory flows were observed at two to three 
times the fill time well before the reflection process could have an effect. 

Future studies interested in precisely defining shelter flow should 
investigate the circulatory flow phenomena. The investigation should include 
the dynamics of the flow process and the energy and momentum transition and 
dissipation processes. From these investigations, estimates of the damage 
potential of these flows should be possible. 

GENERAL COMMENTS 

In reviewing the literature, analyzing the numerical and experimental 
data, and in some of the preliminary analyses, some interesting observations 
on the shelter flow field were encountered. The first of these grew out of 
observations of early numerical simulations. As the percent opening (that is, 
the area of the entrance relative to the wall in which it exists) increases 
over 10 percent, a significant reduction in the reservoir pressure occurs. 

This decrease was documented experimentally as discussed in Appendix C and 
also in the numerical simulations as discussed in Section 3. 

A preliminary investigation indicated that this corresponds to the 
difference between the free stream inflow into the reservoir and the high 
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velocity inflow into the chamber. The difference can be shown by simply 
considering a duct with an orifice. A 5 psi overpressure step shock will be 
assumed to propagate down the duct towards the orifice which has a 20 percent 
opening. To illustrate the difference between inflow and outflow, the mass 
rates can be estimated using the tables and graphs in Refs. 2 and 3. 

On striking the orifice, the 5 psi wave would reflect to 11.4 psi over¬ 
pressure. The reflected overpressure forms the reservoir for the inflow 
into the chamber. Under these conditions, the influx of mass into the 
reservoir can only be as great as the mass rate of the incident shock, while 
the efflux of mass into the chamber is a function of the expansion conditions. 
For the 20 percent opening, the ratio of influx to efflux mass rates would 
be 1.40. This ratio states that 0.4 of influent mass rate is available for 
the reservoir which is expanding at the reflected shock wave’s propagation 
velocity (1091 fps). To maintain the reflected condition, the retention 
mass rate must equal 1, hence a decrease in the reservoir pressure must occur, 
reducing the influx into the chamber. 

This result explains why the chamber-filling computer programs have an 
increasing error as the percent opening increases. The inclusion of the 
work term in the equations adjusts for the work that must be done on the gas 
internal to the chamber, but an allowance is not made for the effect of the 
inflow on the reservoir pressure. Similar results were noted for the trans¬ 
mission of shock waves into tunnels (Ref. 42). As the reflecting surface 
against which the incident shock wave impinged was reduced, the driven shock 
wave overpressure in the tunnel was reduced. 

In terms of shelter flows, this result would indicate that the standard 
clearing time equation overestimates the driving pressure. The rarefaction 
waves generated by the inflow process contribute to the reduction in reflected 
pressure, hence a reservoir pressure value would most likely be between 
reflected and stagnation pressure in the early portion of the filling process. 
The later portions of the process would be "filled" from stagnation conditions 
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and perhaps a lower pressure depending on the mass rate flows into and through 
the shelter and shelter building. 

Further work should be done in this area as the indications are that 
present techniques may overestimate the reservoir pressure. Since the reser¬ 
voir is the driving force for the inflow, lower inflow velocities may actually 
occur. Based on the limited evaluations undertaken in this study, a value 
20 to 50 percent of the difference between stagnation and reflected pressures 
greater than stagnation pressure would appear to be reasonable for estimating 
the reservoir pressure. The future work should relate the phenomena to the 
distribution of open area in buildings so as to be able to provide more 
definitive statements. 

In the review of the literature, reference to swirl of flows in ducts 
after and in turns in the duct were found (Refs. 59, 60, 61). Turns in the 
entrances to shelters are to be expected, for example, stairwells. If the 
inflow does swirl, a stronger lift component may exist than for direct axial 
flows into the shelters. Although a low priority consideration, future 
work might ascertain the possibility or the extent of this type of flow. 

A review of the recent OCD work units on structural failure and debris 
formation provided a different view on the potential for blast transmission 
through elevator shafts into shelters. The results of these studies indicated 
that the period for an elevator door to fail and move a significant distance 
might approach the filling time of the chamber. If this can be substantiated, 
then the filling of the shelter would be solely through the open entrances. 

The above comments do not include the situation in which the elevator may be 
between the aboveground floors and the shelter. 

The last observation to be made on the shelter flows is the temperature 
of the gas within the chamber. An adiabatic expansion into the room and 
adiabatic conditions for the flow in the room are generally assumed. Under 
these conditions the temperature for the room after filling would then 
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correspond to an adiabatic compression to that pressure or approximately the 
temperature behind the incident shock wave. For a 15 psi overpressure shock 
wave, this temperature would be approximately 200°F. In terms of the long 
durations anticipated in the prototype situation, this temperature may be 
sufficiently high when coupled with the air movement in the shelter to 
warrant further consideration. 

Some documentation of the occurrence of this potential hazard can be 
found in the summarization of "Non-Line-of-Site Thermal Problems" in Ref. 62. 
Temperature measurements taken in shelters away from the shelter entrances 
during the 1953 and 1955 test series at the Nevada Test Site did show signi¬ 
ficant temperature increases. Many animals exposed in these tests (including 
the 1957 test series) were also reported to have shown skin and fur burns. 
These burns were reported to have been intensified on the surfaces of the 
animals facing the inflow. As pointed out by Dr. White, these tests were not 
conclusive, but evidence does exist to warrent a further review of the hazard. 
The potential hazard in slanted shelters would be intensified as the expected 
durations would be considerably longer than the duration of the tests from 
which the data was collected. 
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Section 5 

REVIEW OF NUMERICAL SIMULATIONS 

A series of numerical simulations was prepared as a part of this work 
unit at the OCD computer facility at Olney, Maryland, and the computed results 
stored on magnetic tape. Selected parameters were printed out for review, 
and a yet smaller number were plotted in spatial contours for about four time 
intervals per computer run. Even with the selected output, the volume of 
data could not be presented within the economic constraints of this study. 

The purpose of the simulations was to provide further insight into the 
flow process within a shelter to aid in the development of the conceptual 
model. Though the full set of data cannot be presented, a verbal description 
of the results will be provided in this section. The results stress inter¬ 
pretation of the calculated data as it pertains to understanding the flow 
in the shelters. 

CHAMBER LENGTH VARIATION SIMULATIONS 

The first series of simulations undertaken covered shelter spaces with 
10 percent opening on a wall normal to the direction of propagation of the 
incident shock wave. The entries were located at one end of the upstream 
wall. The geometry of the three rooms involved was 30 x 30, 30 x 60, and 
30 x 180 (width x length). The incident shock wave had a 4.9 psi over¬ 
pressure in all three cases. 

The expansion of the shock wave into and through the room generally 
followed the pressure distribution model described in Refs. 2 and 3, A 
strong decay in pressure across the chamber was observed behind the shock wave 
as the shock wave reached the rear wall of the chamber. At this time a low 
pressure area of zero psi overpressure was observed at about five entry widths 
from the front wall. The pressure contours at this time were approximately 
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perpendicular to the side walls of the chambers. This coincides with an observa¬ 
tion that the shock wave first expands, and then is reinforced by the inflow into 
the chamber and retains its approximately normal position to the side walls. 

The pressure distribution throughout the chamber varied spatially with 
time. Local high and low pressure areas formed and propagated through the 
chamber during the filling process. The chamber appeared to fill to about 
75 percent of the exterior pressure in two transit times of a sound wave 
across the chamber (one filling cycle as defined in Ref. 2) . The external 
or reservoir pressure was below reflected pressure (11.2 psi) and approached 
a steady value in the range of 9.5 psi. This result supports the argument 
for reduced reservoir pressures presented in Section 4. 

The dynamic pressure contours for the chambers lagged behind the pressure 
contours. The dynamic pressure profiles were just beginning to approach 
their maximum intensities with their leading edge about 5 to 10 entrance 
diameters into the room as the expanding shock wave struck the rear wall. 

The profile continued to expand into the room as did the point of maximum 
intensity (about 7.5 psi). The dynamic pressure contours remained along the 
wall on the axis of the doorway throughout the entire flow period. Upon 
striking the rear wall, the flow turned and moved along it until it struck 
the next corner. 

The flow, as it turned the corner, increased the local pressures corre¬ 
sponding to a stagnation zone. These local stagnation pressure zones and 
the flows causing them extended well past the filling time and were observed 
to continue past the point at which the pressure in the room had equalized. 

The flow intensity, as monitored by the dynamic pressures, was reduced after 
it turned each corner. The dynamic pressure contours were identical for each 
of the modeled rooms up to the point the rear wall was sensed, 

ENTRANCE LOCATION EFFECTS 

The next series of simulations was comprised of single rooms but their 
entrances were positioned other than in the center of the front or reflecting 
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surface, A free flow field was provided around the chamber. The chambers 
used in the simulations all were 40 ft long in the direction of incident 
flow and 20 ft wide in the transverse direction. All openings on the front 
and rear walls were 2 ft and on the side walls 4 ft. These openings were 
positioned in different combinations in each of the simulations. 

The first simulation had the opening positioned on the side of the 
shelter at the upstream corner. The inflow into the chamber was initially 
similar to that observed in models filling from the front, but the flow was 
skewed toward the rear. The skewing of the flow corresponded to the component 
of the external flow past the entry. Both the flow and pressure distribution 
throughout the chamber were centered on the diagonal from the entry to the 
opposite corner. 

Circulatory flows were not observed in any significant magnitudes. The 
reason for the absence of the circulatory flow correlated with the observation 
that the peak fill pressure monitored in the simulations was 1.5 psi, well 
below the incident shock wave f s overpressure. The entry, located at the 
corner of the shelter, was below the separated free stream line in which area 
a low pressure zone exists. The pressure monitored in this region corresponded 
to the filled pressure. 

The second simulation in this series had the entrance also located on 
the side wall, but at its mid-point. The chamber pressure at the "filled” 
condition was 4.0 psi, which corresponded to the side-on pressure at the 
exterior of the entrance. The maximum dynamic pressure monitored in the 
inflow was 1.5 psi. At the end of the simulation period the inflow had not 
yet reached the opposing wall of the chamber. 

A chamber whose entry was on the rear wall was also simulated but prob¬ 
lems in the plotting made it difficult to use the data. Since the limited 
results gained showed a similarity to the other fill processes reported in 
Ref. 2 and also in this section the simulation was not replotted. A fill 
pressure in the chamber of 5 psi was monitored, which corresponds to the 
free stream side-on pressure. 
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The next chamber simulated in this series had entrances both at the 
front and rear of the chamber. The main inflow occurred through the front 
entrance. Inflow from the rear opening was initiated, but it quickly 
subsided. The pressure in the chamber peaked at the 2.5 to 3.0 psi region 
as a through flow was evolved in the model. The pressure on the front of the 
model was approximately 5 psi. A peak dynamic pressure of 6 psi was monitored 
which steadied to a peak value in the through flow of 5 psi. These intensities 
are peak values and do not represent the average dynamic pressure of the flow. 
Due to the complex geometry of the dynamic pressure contours, good average 
values could not be determined. 

The last configuration tested in this series had openings on the front 
wall and in the center of the side wall of the shelter. The inflow process 
through both openings was initiated separately. The flow in the front was 
accordingly much stronger as the reservoir pressure was greater. The inflow 
from the side was initially skewed to the rear of the chamber, but as the 
inflow from the front was sensed the side inflow began to skew towards the 
front of the chamber. 

The side inflow did not contribute very much to the filling process as 
is evidenced by the peak dynamic pressures for the two flows — 6.5 psi early 
in the inflow from the front steadying to 4.5 at a later time, and a 0.5 psi 
maximum for the side inflow. The pressures in the chamber ranged between 
4 and 6 psi overpressure which settled to a filled overpressure of 4.5 psi. 

This internal pressure corresponded to a pressure on the front of the chamber 
of 5 psi. Similar to other chambers with front entries, very low pressures 
were monitored during the early inflow process at a point just inside the 
front entrance along the front wall (a negative 1 psi was the lowest value 
shown in the plots). 

MULTIPLE CHAMBER SIMULATIONS 

The last series of simulations was undertaken for a two room configuration. 
The chamber was 30 by 90 ft with 3-ft openings positioned on the exterior and 
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the divider walls. The simulations were divided into two sets, the first 
set had the axis of the doors aligned and the upstream inflow occurring in 
the same cross-section as the chamber; the second set had the doors offset 
and the exterior of the chamber extended for 30 ft on either side of the 
entrance to approximate an infinite reflecting surface. 

The simulations of the aligned-door set showed a pronounced through flow. 
The upstream chamber of the first simulation was 30-ft long and the downstream 
chamber was 60-ft long. The upstream room did not appear to fill quickly. 

The exterior pressure reached a steady value of 10 psi. 

The second simulation of this set divided the chamber into two equal 
45-ft rooms. The early flow and pressure histories observed corresponded to 
those of single room chamber. As the filling process continued, the pressure 
in the front chamber filled to 2.5 psi and leveled off while the flow contin¬ 
ued on into the second room. Very little of the inflow appeared to remain in 
the first room. As the second room filled to about 4.5 psi, more of the 
inflow remained in the first room. The pressure in the second room increased 
to a peak of around 10 psi and then dropped back to about 8 psi. 

The peak dynamic pressure monitored in the flow in both rooms ranged 
between 6.5 and 7.0 psi. Circulatory flow patterns were observed to occur 
in both rooms. The circulatory flows were more severe in the second room. 

The strength of the through flow can be envisioned by considering the stagna¬ 
tion regions that formed in the corners of the rooms with solid rear walls. 

In the case of the shelter with an entry to a second room at this location, 
the momentum is transferred into the second room rather than forming the 
stagnation zone. 

In the second set of simulations, those with the offset doors, the same 
two floor plans were tested. The front room in the simulation of the 
30 x 60 ft combination filled quickly to 3.5 to 4.0 psi. The second room 
began to fill from this pressure differential as the inflow from the front 
entry curled around the corner and moved along the partition wall. A 
well-developed flow is formed along the side and partition wall which is 
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largely shunted into the second room. Only a small portion of this flow 

remains in the front room to establish the circulatory flow pattern. The ■ 

second room fills to its peak pressure at a much later time. At 200 msec 

(the end of the computation time) the rear room overpressure was 3.0 psi, the 

front chamber had an overpressure of 5.5 psi, and the exterior a 10.5 psi 

overpressure. 

I 

The last simulation undertaken in this series had two 45-ft chambers 
with offset doors. The flow pattern appeared to be similar to that for the 

! 

30 x 60 ft condition. The pressure in the front chamber increased to ] 

4 to 4.5 psi, as the pressures in the second chamber rose to 2.5 psi at 

i 

about midway in the computation of the fill period. The computation was j 

ended at this point due to the long running time of this particular type of 
problem. The dynamic pressure monitored in the inflow in the second room 
had a peak value of 3 to 3.5 psi. The outside pressure again leveled off 
at 10.5 psi overpressure. The pressure monitored in the corner of the first 
chamber, near the offset entrance, was slightly reduced due to the flow into 
the second room, particularly during the early portion of the flow process. 

COMMENTS 

The preceding review described the relationship between the driving 
pressure and the peak dynamic pressure. The relationship developed in Ref. 2 
is presented in Fig. 39 with the data points for peak values outlined in the 
discussion plotted over the curves. These peak conditions appear to follow 
the free flow vena contracta graph and generally fall near the upper portion 
of the bounded region developed in Ref. 2. 

The values discussed are peak values and do not represent the space-time 
history of the inflow process. These representations are quite complex and 
at the conclusion of this work unit were not in suitable form for presentation. 

It is planned that a further attempt at a representation of the flow will be 
made in the companion work unit to this study (1154G). 

1 
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Fig. 39. Average Entry Dynamic Pressure 
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Section 6 
DISCUSSION 

One objective of this study was to develop a conceptual model of the air- 
blast-induced flow into basement shelters. The term conceptual model was 
not defined in the scope, other than to indicate it should be useful in 
defining engineering estimates of the flow. Before further considering the 
flow, some implications of the term conceptual model will be reviewed. 

The first step in reviewing the term ’’conceptual model” was to graphically 
portray the flow field boundary conditions. The flow chart developed in 
Fig. 40 is often referred to as a system diagram, hence the process portrayed 
will be viewed as a system using the standard nomenclature of general systems 
theory. The purpose of this exercise was to graphically portray and explore 
the relational structure between the components of the flow-shelter inter¬ 
action upon which any model must be built. 

The system will be defined as the flow field into and within the chamber, 
its subsystems being the inflow and the flow internal to the chamber. The 
environment for the system is made up of the shelter area and reservoir 
conditions. The suprasystem includes the surroundings, the air blast environ¬ 
ment, and the building in which the shelter is located. 

If these designations are reviewed, the elements of the system components 
can be determined. For example, the shelter area (part of the environment) 
has two elements — entries and interior. These elements in turn have 
properties which were defined as location, size, and geometry. A hierachical 
structure can thus be defined. The simplest prediction procedure would result 
if each of the levels in the structure hierachy were independent and could be 
considered individually. 

The condition of separable components is equivalent to the process of 
superposition of parts in a linear system. For these conditions to exist, 
each component of the overall system must be independent. An alternative 
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Fig. 40. A Simple Graphic Representation of Flow Field Boundary Conditions 
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view of the problem which provides some insight into the process would be 

to consider the systems diagram as a set S comprised of a group of subsets 

S , S , ., S . If an element e of .is also an element of S_ , then 

1 2 n 1 2 

S and S have an intersection. The effect on element ^ by the properties 
1 2 

of S must then affect set S , as element e is common to both. If this 
f 

interaction occurs there is a relational structure between sets and they 
cannot be considered independently. 

Two particularly strong interrelationships exist in the flow process in 
shelters which can be seen as feedback loops in the system diagram. These 
feedback loops occur between the inflow and the reservoir conditions and 
between the flow in the shelter and the inflow. The subsystems, as defined 
above, are therefore interrelated, as is the environment to the subsystems. 
These processes cannot, therefore, be considered to be independent with 
solutions that can be superpositioned on each to form a general solution. 

An effective conceptual model must therefore not only define the 
individual subsets and environment components, but also their relational 
structure on the basis of a continuing interaction. To define the inflow 
into the shelter the details of the evolving internal flow fields must be 
established. These very technical concepts must then be interwoven to form 
the description of the total flow field. The complexity of this weaving pro¬ 
cess was stressed in Section 4, but more importantly it stresses the overall 
problem faced—i.e., too many potential variables for any explicit solution to 
be developed. 

For the purpose of argument the flow field will be assumed to be defin¬ 
able in space—time—intensity parameters. The problem then simply reduces to 
defining the suprasystem and environment elements and their properties. The 
shelter area has a direct input into the flow system, therefore the definition 
should be precise if the conceptual model is to be meaningful. For example, 
two questions that must be answered are: 1) What specifically is a blast- 
slanted shelter? and 2) What properties does the shelter exhibit? Without 
specific answers to these questions it is difficult to describe a model of 
the flow, as the environment for the system has not been defined. For an open 
system this lack of definition means an undefined system exists. 
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The properties of the environment and suprasystem must also be defined. 

A preliminary review of previous work on the open shelter concept and of 
the RTI survey of new construction (Ref. 63) showed that a wide distribution 
in shelter geometry, size, and location within a building occurred. Based 
on this survey a blast-slanted shelter cannot be described explicitly, but 

/ rather can only be defined in terms of statistical distributions. Work on 

/ 

evolving more complete bounds on the definition process was begun, but time 
constraints precluded its continuation. The preliminary survey material was 
transferred to Work Unit 1154G in which a thorough analysis of the shelter 
data of the RTI survey is being made. 

Another factor which stresses the problem of interaction between 
components is the orientation of the structure and the blast wave and the 
entry location in the interaction zone. In Sections 3, 4, and 5 observations 
and conclusions were discussed which stressed the strong role these inter¬ 
actions play on the inflow and internal flow processes. In developing a 
conceptual model the possible variations in the interaction zone must be 
defined, i.e., building—biast wave orientation, blast wave properties, 
building geometry, shelter location, and entry location. In the earlier 
review of the system it was noted that the inflow affects the reservoir 
conditions as do the above; hence an explicit definition again cannot be 
undertaken. 

Rather than belabor the difficulties associated with the definition 
process, the positive aspects can be considered. Based on the preceding 
arguments a model can be defined, but only in general terms. These terms can 
verbally describe the processes, but do not, as such, provide a relational 
structure from which quantitative estimation procedures can be developed or 
applied. 

A general model of the inflow can then be described only for a specific 
situation. The situation most commonly selected is for the incident shock 
wave to reflect normally from the shelter wall containing the entrance to 
the shelter. The process for the air-blast-induced filling of a chamber was 
discussed in Refs. 3 and 4 and then reviewed in Section 4. In Section 4 
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shortcomings in the early descriptions were noted and discussed. Rather 
than reproduce the material describing the inflow and the internal flow 
processes the summary will stress how these descriptions, when combined with 
the discussions of Sections 3 and 4 and Appendix C, can be used to provide 
an estimate of flows in the shelters, and to further illustrate the difficulty 
of developing a general conceptual model. 

The first approach to defining the flows in the simple chamber is to 
define the reservoir conditions. The gross reservoir conditions are defined 
by the blast wave structure interaction. To refine the accuracy of these 
conditions the inflow into the shelter and into the structure must be evolved. 
The effect of the inflow can then be computed for each specific situation. 

The predictions can be made either by analytical, experimental, or 
numerical procedures. The most likely procedure to be used would be the 
analytical approach. Using as a guide the requirement of engineering accuracy 
in predictions the following procedure is suggested. Assume the building in 
which the structure is situated is a solid block and calculate the reservoir 
pressure history for the shelter using standard procedures. If the openings 
on the upstream surface comprise less than 10 percent of that surface area 
use the calculated values. 

If the openings on the building’s exterior exceed 10 percent, the 
reservoir pressure should be adjusted. It is suggested that the adjustment 
use the methods outlined in Ref. 64. The following formula (which, based 
on observation, assumes a square power distribution) provides an approxima¬ 
tion of decreases in reservoir pressure for long duration waves: 

P a (t) = p(t) + U - ) [ P(t) " P s (t) ] (32) 

where p(t) is the calculated reservoir pressure, p g (t) is the fr ee stream 

stagnation pressure, p (t) is the adjusted reservoir pressure, and PO is 

a 

the percent opening on the front surface of the building. The preceding 
assumes an equal distribution of openings on the surface. Even if this 
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condition is met, the above equation is an approximation, but it should 
provide estimates within a ± 10 percent variation. 

Using the calculated reservoir conditions the inflow into the chamber 
can be determined from the graphs in Ref. 2. These graphs will provide the 
maximum inflow. An acceleration of the inflow from the flow velocity behind 
the incident shock to the peak velocity will occur in approximately five 
sound wave traverse times across the entrance. The location of the peak 
velocity will range between 2 to 6 diameters from the entry, depending on 
when in the inflow process it has occurred. 

The peak velocity can be assumed to be constant for the fill time of the 
chamber. For estimation purposes the peak velocity can then be assumed to 
decay linearly to zero in another fill time. The preceding provides a guide 
for establishing the peak intensity history of the flow, but it does not 
provide the spatial intensity history. The difficulty of quantitatively 
modeling the flow in space—time terms, not only in the transitory stages but 
also in the long term flow patterns, was brought out in earlier discussions. 

Based on recent work the basic character of the flow can be established. 
The intensity along the inflow axis decays exponentially with distance between 
the location of peak intensity and the rear wall. The location of the peak 
intensity moves inward with time. The increase or decay in intensity with 
time at a given location varies with the location. The intensity of the 
inflow also varies in a transverse direction. The combined effects make it 
difficult to undertake to predict the space—time—intensity of the full flow 
field. 

The prediction of the turning process as the flow strikes the rear wall 
is dependent upon the ability to predict the conditions of the impinging jet. 
From the numerical simulations undertaken the dynamic pressure of the turned 
fluid appears to be one-half of the incident intensity. Part of the reduction 
was due to the accumulation of fluid in the forming stagnation zone and, in 
later times, to the general decay of the inflow process. 
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The above discussion was for a single chamber, with a simple geometry, 
without multiple openings, and for constant conditions in the reservoir. 
Inclusion of these variables make the description of the flow yet more com¬ 
plicated. These comments support the earlier conclusion as to the difficulty , 
of establishing a generally valid conceptual model of the inflow. The in¬ 
ability to document the vast majority of information and observation on the 
flow process developed in this study, due to space limitations, further 
stresses the lack of an adequate description of the flow. 

A general understanding of the overall scheme of biast—induced flow 
processes within blast-slanted shelters exists. An adequate though general 
understanding of the individual flow processes also exists, with some further 
study required in only a few specific technical areas. Based on the results 
of this study it appeared that further studies of technical areas should be 
in direct support to the specific needs of some facet of the shelter research 
program. The work undertaken in this study also indicated that a thorough 
all-inclusive and general conceptual model from which engineering estimates 
may be made for any situation is not a practical goal. To produce a model that 
would simulate flow fields and pressure distribution for any geometry with a 
90 percent accuracy for the full blast profile would require inordinately 
high expenditures and amounts of technical effort. 

The needs of the shelter research program do not appear to require this 
level of technical understanding. The flow as described, with some improve¬ 
ments in the understanding of a few specific technical areas, should be 
sufficiently accurate to answer specific questions. The specific regions of 
flow that do require additional understanding are the long term history of 
the jet in the chamber, and how it forms into the circulatory flow. The 
dissipative mechanisms and energy redistribution within the shelter should be 
included as part of this investigation. 

Another flow area that should be investigated is the effect of the inflow 
into the shelter, and the interior of the shelter building, with respect to 
its effects on the reservoir pressure. Closely allied to the flow is the 
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effect of the building shock wave orientation, the geometry of the building, 
and the location of the shelter. In terms of the system diagram discussed 
earlier, the problem again translates to the definition of the environment 
components which form the inputs for the inflow processes. 

Further arguments in support of considering these system constraints 
in more depth can be evolved from the descriptions of the flow in Section 5. 
The simulation discussed in Section 5 had the openings located at different 
points around the periphery of the shelter. The alterations in geometry 
were shown to drastically affect both the filled pressure and the peak inflow 
intensities. These effects could also be related to the reservoir conditions. 
The monitored velocities and pressures in these simulations were directly 
related to the reservoir conditions and were predictable from present adia¬ 
batic expansion theory once the reservoir conditions were defined. 

The side and rear entrances of the multiple entrance configurations 
were seen to add very little to the inflow into the chamber. The rear 
entrance reduced the air circulation within the room as it formed a through 
flow across the shelter. The results of the simulation do indicate that the 
model of the inflow as described for the simple shelter does break down. The 
results described were for one angle of shock wave incidence relative to the 
shelter. Due to realities of construction and the inability to relate 
building orientation with the detonation point, the reservoir conditions 
cannot be defined — making a general inflow model harder to define. 

The multiple room numerical simulations showed a strong dependence on 
internal geometry as well as reservoir conditions. The location of the 
divider wall determined the circulatory flow patterns, how quickly they 
form, and to what pressure level the first room fills. In the situations 
with aligned entries the shelter geometry with the small upstream room showed 
the strength of the through flow was sufficient to carry all the inflow 
through the first room. This situation was similar to the front-rear exit 
simulation of the single room. 
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The multiple room experiments again showed (Refs. 2, 3, 4, 5) that the 
upstream room does not completely act as a reservoir for the flow into the 
second room. In some cases the flow was seen to move completely through the 
upstream room and not add appreciable mass to the first room. The fill in 
most cases appears to be through the upstream room to the rear chamber and 
then back to the upstream chamber. Based on this model the downstream room 
will always fill to a higher pressure. Situations that were not investigated 
were; small downstream rooms, multiple openings to rooms, and varied incident 
shock wave to shelter orientations. 

A question can be raised from the above results; why does the upstream 
room pressure increase if the flow is through the upstream room. The explana¬ 
tion appears to be related to a combination of the piston effect, discussed 
in Section 4, and the dynamics of the shock wave propagation through the 
room. The shock wave, as it expands through the room, moves gas from the 
vicinity of the entrance corridor. In a normal situation of an expanding 
shock wave, the rarefaction wave decreases the pressure behind the shock wave 
and the negative trough of the wave initiates a flow in the opposite direction. 
The reverse flow would be supported by reflected waves propagating back 
through the chamber. 

The inflow, however, occurs behind expansion processes. The inflow, as 
discussed earlier, drives the shock wave and hence spreads the shock front 
velocity, or alternately, it strengthens the expanding shock wave. The 
reverse flow cannot, therefore, flow back into the area in and about the 
inflow corridor from which it flowed originally as the inflow is now occupying 
that space. It is felt that this simple process explains the early rises in 
pressure in the chamber before the influent mass can raise the pressure 
throughout the chamber. 

The most important aspect of the flow review undertaken in the work 
unit was that it pointed some ways of planning future work in support of the 
OCD shelter research program. The work stressed a need for a more systematic 
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approach to the research. A thorough understanding of the flow can be 
gained with added effort, but the depth of this understanding can not be 
determined without a stronger emphasis on the practical aspects of the 
problems faced. Using the practical considerations as guides it may be 
possible to avoid some of the difficult technical questions. 

The framework of understanding of the flow processes is in general 
sufficiently precise to predict flows with engineering accuracy if there 
is some idea of what the shelters might look like and what external condi¬ 
tions can be assumed. Assumptions of the worst conditions, although con¬ 
servative, establishes a rather poor estimate of actual conditions. In 
the post-attack research area, the relational problem has begun to be con¬ 
sidered (Refs. 65 and 66). Also, Ref. 67 stressed the importance of the 
geometric parameters in determining both flows and loading. The study under¬ 
taken by this work unit also indicated that a similar approach to the shelter 
problem might be taken in the immediate future. Under these guidelines, the 
term conceptual model would be used to include all the aspects of the air 
blast shelter interaction, not simply the flow processes. 

From discussions held with other OCD engineers working in the shelter 
research area the above review might also include their specific areas of in¬ 
terest in the shelter research program. Different levels of understanding of 
the flow exist for each of these, but at present a good working outline or 
systematic representation is not available. The three working areas that 
appear to exist are: 

• The pressure loading on structural elements 

• The injury potential in the below-grade slanted basement shelter 

• The above two problems with respect to above-grade building areas 

Of the above, the evolvement of loading profiles on structural elements is a 
simpler procedure than providing general predictions for flow. 

If the above reasoning is taken a step further, the question is raised 
as to whether the flow field is the item of interest, or if the damage 
potential is of interest. From the results of the investigation undertaken 
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in this work unit the definition of the damage potential appeared to be more 
influenced by the external parameters than by improved prediction techniques. 
The purpose of future studies may be to ascertain the damage potential and 
then suggest ways of limiting its causes. Based on the work undertaken a 
simple prediction procedure for providing flows for all possible geometries 
is not feasible, if the tenet of simplicity and small computer size is 
maintained. 

The work on the code in Section 3 showed that for a simple two-dimen¬ 
sional geometry it might be possible to write a simple code. Three-dimen¬ 
sional effects and complex geometries would have to be excluded. The effort 
to provide this form of code would vary with the complexity of the output 
desired and the number of output variables. For the prediction of pressures 
in shelters it might, however, be an economic feasibility. 

The use of the above code, if developed, would appear to be an over¬ 
sophistication relative to the real-world environment. The aim of this work 
unit was to consider the concept of the flows, and it brought out some of the 
oversophistication of approaches that have been used to generate data and 
develop theories. The reality of the prototype situation produces variables 
and stochastic distributions that far outweigh many of the minor fluctuations 
observed and then researched. This problem is not unique to the aerodynamic 
portion of the shelter research program. A pronounced tendency to get involved 
in technical detail , rather than the broad general overview was observed (for 
example, a review of this report will underscore this tendency). 

Based on the work undertaken in this work unit the following suggestions 
are made for future studies: 

• The first priority of future studies in this area should be to 
develop a system for describing an overview of the processes involved. 

• Further studies should work towards the final objectives of the 
shelter research program area using the system model as a guide. 

The use of this approach should significantly increase the 
effectiveness of further studies. 
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• Experiment should be used to corroborate analytical work, the 
systems model, or to provide answers to specific problems. Rather 
than increased experimental effort, more effort should be placed 
into interpreting the results to provide empirical insights and 
knowledge of the air blast shelter interaction. 

• A freer use of analytical and numerical simulations should be 
undertaken to aid in the evolvement of understanding, but under 
the same constraints as with the experimental simulations. More 
emphasis should be placed on full-scale simulation by both of the 
simulation techniques. 

The following specific recommendations on technical areas for further study 
of the flow are made: 

• A brief review of the arguments on the spread of turbulence and 
the growth of boundary layers via different procedures should be 
undertaken to ascertain if the implications derived in this work 
unit will bear the weight of further scrutiny. 

• A further expansion of the piston effect analysis should be made, 
particularly relating the spread of the inflow and the definition 
of the inflow movements and distribution of mass and energy. 

• A better understanding of the transition of the inflow energy 
and its distribution, in conjunction with the evolvement of 
the circulatory flows should be developed. A combination of 
experimental and numerical simulations appears to be needed to 
develop a better intuitive insight so as to shorten further 
analysis. 

• The problems of interactions between internal rooms within 
shelters should be studied. 

• A better understanding of the effect on the reservoir conditions 
of the flow into the shelter or shelter building should be 
developed. 


A closing suggestion is that the work in this area be closely tied into 
operational concepts as defined in Ref. 68 to ascertain what effect these 
will have on the work area being studied. The system model discussed 
earlier would fit into this framework. 
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Appendix A 

REPRESENTATIVE FINITE-DIFFERENCE EQUATIONS 

The following finite-difference equations developed by Rich are presented 
in this appendix (as Rich developed them in Ref. 13) to provide the reader a 
representative set of Eulerian finite-difference equations and calculation 
procedures. In this presentation, the center of each cell, (i,j) is labeled 
by its row (i) and column (j) location in the mesh. Time is differenced in 
multiples (n) of constant time increment (At). The values of density 
velocity fuf 1 ^ , vf H ^l , and specific energy |E f n ^ | are constant within each 

L i,J i » J J L i > J J 

cell for each time increment n as defined by the superscript (n). 


The equations will be given in the computational sequence described in 
the text. This computational sequence had three phases: (I) the calculation 

of intermediate values of velocity and specific energy in each cell with the 
boundary condition of no mass flow out of the cell; (II) mass transport across 
the boundaries based on the intermediate values; and (III) the computation of 
the resultant conditions in the cell. 

The difference equation for the conservation of momentum and energy during 
Phase I is given by: 
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Where the pressure and velocity at the cell boundaries are given by linear 
interpolation as in the following examples, 
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and the time derivatives are defined by 
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Substitution of the above equations in the original difference equations 
establishes the following equations for the intermediate velocities and 
specific energy (denoted by the tilde) in any cell (i,j): 
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The tilde values represent the change in momentum and energy due to the 
surface forces acting on the control volume, i.e., any cell (i,j). These 
values are used to compute the mass transport effects in the following Phase II 
equations. The mass transport across a cell boundary is proportional to the 
fluid density and the normal component of velocity at the interface. The 
exchange must conserve mass, hence the Eulerian expression for this condition 
must be satisfied and is expressed in the difference format as 
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AxAy - 
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Similarly for the cell boundaries (i-l/2,j) and (i,j-l/2). 
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In the above expressions, the cell boundaries across which the mass is 
exchanged are indicated by the subscripts. 


The tilde values in the preceding equation are those at the cell bound¬ 
aries , not at the center of the cell. Rich (Ref. 13) and other investigators 
have shown the use of simple averages for the densities and velocities at the 
cell boundaries leads to minor, but noticeable, instabilities in the solution 
of the difference equations. Rich's technique avoids the instability by 
preferentially weighting the density and velocity in favor of the cell out 
of which the fluid flows, and to take advantage of the differencing property 
which limits the center of the fluid element (located at the center of the 
cell of the Eulerian mesh) from leaving the cell during the computational time 
step. A better representation for the values of the boundary quantities 
(u, v, p), particularly in the mass flow calculations, were shown to be 
obtained from the first two terms of a Taylor expansion about the center of 
the emptying cell. The use of simple averages for the density terms was shown 
by Rich (Ref. 13) to be sufficient in most instances. 


The Taylor series expansion is about the center of the emptying cell, 
thus the equations for the boundary quantities are dependent upon the direc¬ 
tion of flow at the cell interface. The selection of direction can be best 
accomplished by using the quantity [uf n ^ + u^ ] — in essence the average 

L i,j 1+1,jJ 

value — as the test condition for direction. The average insures the same 
condition is used as a test irrespective of the direction at which the compu¬ 
tation is approached, i.e., either from i or i+1 the result is the same. The 
mass flow ' would be calculated for the right-hand boundary of cell 

(i,j), i.e., (i + 2 , j) , as follows: 
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The flows across the other three sides of cell (i,j) are obtained in a similar 
manner. 

In the Eulerian calculation scheme, the mass flowing from one cell to 
another carries with it the tilde velocity and specific energy of the cell 
from which it came. At the completion of the mass transport of Phase II, the 
composition of a given cell may be represented as a summation of all the mass 
and energy flows across its boundaries. The determination of final values is 
the third phase of the computational sequence. By combining the laws of con¬ 
servation of momentum and energy with the mass transport quantities, the ter¬ 
minal value of velocity and specific energy for each cell can be obtained: 
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To simplify the computational scheme, Rich numbered the four sides of 
cell (i,j) k = 1, 2, 3, and 4, as shown below. 


4 



A function C. . (k) for any cell (i,j) and its side k was defined such that 
i > J 

C. .(k) = 1 if fluid flows into cell (i,j) across side k 

1 7 J 

= 0 if fluid flows out of cell (i,j) across side k 

Rich then presented the following difference equations for the final velocity 
and energy, 
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(Equation continued) 
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These equations complete Rich's computational cycle for the transition of the 
fluid configuration from time n to time (n+1). 
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Appendix B 

A DISCUSSION OF THE GEOMETRIC ASPECTS OF PRESSURE MEASUREMENT 

The wide use of pressure—time measurements from scaled experiments in 
the Office of Civil Defense (OCD) shelter program requires some knowledge of 
the characteristics of the mapped pressure histories. The following discus¬ 
sion considers some aspects of the pressure-mapping process which have not 
been widely discussed in other sources. The material is based on an unpub¬ 
lished draft technical note the author prepared while employed at the Ballistic 
Research Laboratories (BRL). The gages and data referred to in the discussion 
are those of the BRL. Since the main body of available data on flow within 
OCD shelters has resulted from scale model studies done by the BRL, it seemed 
in the best interest of the OCD shelter program to leave it in that format. 


INTRODUCTION 

The aerodynamic effects on a prototype object that are created by its 
interaction with an air blast wave can be experimentally determined by testing 
linearly scaled models in shock tubes. One measureable aerodynamic property 
of the interaction field surrounding the object is its overpressure as a 
function of time. To obtain pressure histories related to points within the 
field, pressure transducers are mounted on and about the modeled object. 

A number of these localized histories can be combined to establish the over¬ 
all aerodynamic loading on and about the object. 

The developed loading on the model can be related to the prototype situa¬ 
tion by the use of scaling laws. In this appendix, linear scaling will be 
used, i.e., there will be a one-to-one relationship in aerodynamic properties. 
Flow velocities are then equal, model-to-prototype situation; and, therefore, 
time and distance scale linearly. The scale factor of the system is the ratio 
of the dimensions of the prototype to those of the model. 
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As the scale factor increases, i.e., the dimensions of the model decrease, 
the fixed size of the gage’s sensing element becomes larger in relation to the 
model's surface. In linear scaling, corresponding to the decreased dimensions, 
there is a decrease in the time scale of the aerodynamic events being measured. 
These factors make the gage response, the average of the pressure acting on its 
surface, more the history of an area than of a point. The shift is due to the 
larger relative area of the gage, in comparison to the modeled surface, as well 
as to the shortened time scale of the aerodynamic occurrences. Measurement 
error is introduced due to the fixed velocities having to cross longer scaled 
distances, while the properties of the field vary more rapidly with both dis¬ 
tance and time. 

One such condition occurs with the beginning of the mapping of the pres¬ 
sure upon incidence of the pressure wave at the leading edge of the transducer 
element. As the wave traverses the gage, an increasing area is exposed to the 
pressure until the pressure wave has completely crossed the gage element. At 
any given moment, the pressure mapped by the transducer is the average of the 
pressure then acting on the sensing element. If the wave changes its amplitude 
before crossing the gage, as occurs with sharply decaying exponential waves, 
then the true maximum pressure is never mapped. 

In the following discussion, this form of measurement will be considered 
for the linearly and exponentially decaying shock waves, such as those which 
are formed in the interaction of shock waves with scaled shelter models. The 
interaction forms a complex wave system, which can be considered to consist of 
a series of waves, as shown in Fig. B-l. For this investigation the shortest 
duration of any component wave will be 50 jusec for a linearly decaying wave, 
and 100 fxsec for an exponentially decaying wave. These durations are compatible 
with the test conditions of experimental model shelter studies made at the 
Ballistic Research Laboratories for OCD. 

The decay of the amplitude of the pressure pulse begins at each point on 
the gage at a rate that is measured in time from incidence of the wave at that 
point. For a rapidly decaying wave, the time delay due to the transit time of 
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Fig. B-l. Component Wave Definition 
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the wave over the gage would result in a measureable nonuniformity in pres¬ 
sure. This nonuniformity would be absorbed in the averaging process and would 
result in the following: an inability to measure the wave's peak pressure, a 
phase shift in the mapped profile with respect to time measured from the wave 
incidence on the leading edge of the gage, and an undetermined line of action 
of the average pressure. 

These three areas will be considered within the bounds of the conditions 
that have been set. A mathematical model of the physical situation will be 
established. This model will then be used to show how the mapped pressure—time 
profile is affected by the mapping process. A corrective procedure will then 
be postulated. 

DISCUSSION 

A mathematical model of the gage sensing area's interaction with the pres¬ 
sure wave will be used to represent the mapping of the applied pressure. Let 
a function p(t) be the overpressure-time profile of a shock wave that traverses 
the gage’s circular element of radius R at a velocity U. The function p(t) is 
assumed to vary only in time in the direction of the wave propagation, as 
illustrated by Fig. B-2. 

The distance x is measured across the gage element from the leading edge 
of the sensing surface, and time t is measured from the wave's incidence at 
that edge. The time delay between incidence of the wave front at the leading 
edge and its reaching any point is then x/U. In this coordinate system, the 
pressure at the leading edge begins its decay at t = 0; but for any other point 
x i , the decay begins at ^ = x./U, as illustrated by Fig. B-3. The different 
starting point in time for the decay of pressure (phase shift) is the reason 
the transducer cannot precisely map the incident peak overpressure. The 
valid time range for this mathematical model will begin when the incident wave 
completely crosses the element and will end when the range of p(t) ("the positive 
duration") is terminated at the leading edge. The mathematical model presented 
has the advantage of being spatially and temporally continuous in the defined 
ranges. 
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P(t) 


Fig. B-2. Incident Wave Crossing Gage Element 
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As discussed, the mapped pressure at any time t is the average of the 

pressure acting over the surface of the sensing element. The average, or 

transduced signal, can be described by the integral of the applied pressure 

over the area of the sensing surface. The integral is that of a force volume 

(the depth dimension being pressure), with the property that when the volume 

is divided by the area of the sensing element, the measured average pressure 

is acquired. Figure B-4 presents the force volume as a truncated circular 

* 

cylinder in Cartesian coordinates. Integral calculus (Refs. 1, 2) defines 
the force volume by, 


R x (y) 
r /*2 


V f = 2 


If 

o x (y) 


f (x,y) dx dy 


(B.l) 


where 


f(x,y) = P 


(* - !) 


(B. 2 ) 


x g (y) = R + JR 2 - y‘ 


(B. 3 ) 


x (y) 
1 


R - % IT 


(B. 4 ) 


This integration was carried out for several exponential decaying waves and 
for a linearly decaying wave. A general solution was determined of the form 


p(t) = p(t - R/U) 


f (wR) 


(B. 5 ) 


where wR = cR/U and p(t) is the measured pressure. The term cR/U is the ratio 

3|e 

of one-half of the sensing-area crossing-time to the duration of the wave. 


References for this Appendix are presented at the end of the body of the 
Appendix. 

¥ 

Defined as the time interval 2R/U. 
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The term f(wR) was found to be a series with a first term of unity and the 
second term of order two in wR. The integration for one type of exponentially 
decaying wave is carried out in Exhibit I. In the development of the analysis 
in Exhibit II, it was determined, for the range of experimental values being 
considered, that terms of order two in wR are negligible. With this conclusion, 
Eq. (B.5) can be rewritten as 


p(t) = p(t - R/U) 


(B. 6 ) 


The line of action of the average pressure would be through the centroid 
of the force volume. Exhibit II determines the location of the centroid for 
a triangularly decaying wave through the use of a triple integral. The table 
in Fig. B-5 tabulates the deviation of the centroid from the centerline for a 
ratio of crossing-time to duration of 10. It can be seen that except for 
times past 0.9 of the duration, the deviation is small, with a maximum of 11.1% 
at the end of the integration range of time. 


The results gained in Exhibit II can be generalized by consideration of 
the "force volume 1 ’ as a composite volume. Consider the force volume to consist 
of a cylindrical volume and a wedge-shaped volume (Figs. R-4 or B-6). The 
centroid of the force volume is the resultant of the composite volumes. As 
the ratio of crossing time to duration increases, the portion of the total 
volume occupied by the wedge-shaped volume decreases. With the decrease, the 
effect of its centroid on the centroid of the composite force volume diminishes. 
In the limit, the force volume centroid is that of the cylindrical portion. 


The analysis for the general case by use of centroids of composite 
volumes is carried out in Exhibit III. It is shown that the centroid x of 
the force volume can be given by 



R 



x 
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CR _ 1/2 CROSSING-TIME 
U WAVE DURATION 


1 - Ct 

X 

0.05 

1.111R 

0.10 

1.083R 

0.20 

1.056R 

0.50 

1.028R 

0.80 

1.019R 

0.95 

1.007R 



INCIDENT WAVE GIVEN BY: p(t) = p(l - ct) 


Fig. B-5. Effect of Wave Decay on (x) Centroid 
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where and are the pressures at the lead and aft edges respectively at 
any time t. This result was based on a linear pressure gradient across the 
transducer with the centroid of the wedge at its worst location, 2R. The use 
of the 'worst condition” will provide a maximum value of x irrespective of the 
probable wave shapes to be encountered. 

By this technique an estimate of the location of the centroid may be 

acquired. To use it, one could use the crossing-time and then select a time 

and measure the pressure differential. From this value the ratio p /p can 

2 1 

be acquired and hence the centroid. Within the bounds being considered, the 
decay of pressure does not affect the centroid of the force volume significantly. 
The line of action of the average pressure mapped by the transducer is then the 
geometric center of the sensing area. 

Combining the results of the two analyses, it is found that the pressure 
record is shifted one-half of the gage crossing-time and acts through the 
centerline of the gage. Correspondingly the correct record would be repre¬ 
sented if the record were to be read from one-half the crossing time instead 
of from the first rise of pressure. This record would be valid for that pro¬ 
file sensed at the centerline of the gage. 

The peak pressure is not mapped. This, however, can be inferred by 
extrapolation of the pressure profile back in time to a point equal to one- 
half the crossing-time. If a small flat top wave were included on the top 
of the wave, as shown in Fig. B-7, or if the wave were a compression wave, 
such an extrapolation would be erroneous. The flat top wave is most often 
formed in the regular reflection process or from the compression waves in the 
shock wave formation process. Both these situations are predictable and 
caution can be exercised for these cases. 

Three methods of extrapolation are commonly used: use of the initial 
slope, mathematical curve fitting, and graphical methods. The use of the 
initial slope is treated in Ref. 4, and some inherent mathematical errors 
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Fig. B-7. Pressure—Time Profile of a Regularly Reflected Shock Wave 
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are pointed out. A more pronounced problem is created by distortion in the 
mapped waveform attributable to transducer or the electronic monitoring/ 
recording system characteristics. Since the initial slope method relies on 
tangency at a point, any minor variations can act as levers for producing a 
large error in the extrapolation. 

An example of this is shown in Fig. B-8, The gages were exposed to 
identical waves and were monitored with the same electronic system. The 
differences were the element diameters and the frequency response characteris¬ 
tics of the gage. The BRL-Bar gage is a 2-MHz gage, while the BRL-3 is a 
400-kHz gage. The extrapolated peak pressure varies 6.0%- 

A better extrapolation method would make use of a longer baseline. 

Both of the last two methods described are based on this principle. For 
simplicity's sake, the graphical approach was used in the example shown in 
Fig. B-9. In this figure the same gage records shown in Fig. B-8 were used, 
but the wave forms were extrapolated (by means of a French curve) back in 
time to a point equal to one-half of the crossing-time. The peak pressures 
agree within 1.8%. The dashed curve represents a shift of the upper curve 
to compensate for the difference in time scales, since time was measured from 
the leading edges of the gages. The time shift was determined from the equa¬ 
tion 


R - R 

Time Shift = --—- (B.8) 

or the difference of one-half of the individual crossing times. The remaining 
deviation between the curves is well within the limits of experimental 
accuracy. 

In the preceding discussion, the incident pressure wave was assumed not 
to decay in crossing the sensing element. To ascertain the validity of the 
assumption, an analysis was considered using a wave of constant duration, with 
a variation only in the peak pressure. Based on these conditions, the mapped 
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TIME MEASURED FROM ARRIVAL AT LEADING EDGE OF 
SENSING SURFACE (ysec) 


Fig. B-9. Effect of Sensing Element Size on Mapping Profile 
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pressure profile of a triangularly shaped pressure pulse would be given by 

p(t) = P J 1 - c [t - (i - f)|] j (B.9) 

as developed in Exhibit IV, where (a) is the fractional decay of the peak 
pressure across the gage. For up to a 10% decay in pressure, there is an 
approximate 1% change in R/U, which is not significant to the overall pressure¬ 
mapping accuracy. If a/8 is taken as very small, Eq, (B.9) becomes 

p(t) = p £l - c^t - (B.10) 

This is the equation for the mapped profile of a triangularly decaying wave. 
Based on this result, within the range of interest being discussed, any small 
decay encountered in crossing the element does not significantly affect the 
mapping of the pressure pulse. 

In the discussion to this point, the pressure—time profile has been 
assumed to be symmetrical about the gage diameter in the direction of the 
shock wave propagation. If, at a later time, a second rarefaction wave system 
crosses the gage from another direction and interacts with the initial wave 
system, then that assumption is in error. The line of action of the pressure 
would have to be considered in both the x and y directions. The volume of 
the wedge (Figs. B-4 and B-6) and its effect on the composite x centroid has 
been shown not to affect the line of action of the average pressure. By the 
same method it can be reasoned that the y centroid would not be significantly 
affected by the incidence of a second rarefaction wave system. 

The rate of decay caused by the second rarefaction will not shorten the 
duration of the incident wave. The averaging process described is then not 
appreciably affected. Secondly, the slow rates of decay caused by rarefaction 
systems lend themselves to averaging by the small gage diameters. For these 
reasons it is concluded that in mapping, only insignificant alteration of the 
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pressure profile will be created by the incidence of a second rarefaction 
wave system. If the second wave system were a shock wave, then the technique 
discussed in the earlier portions of this section would have to be reapplied. 

One example of how this analysis can be applied is illustrated by 
Fig. B-10. The finite crossing-time produces the averaging process in the 
mapping with a corresponding error in the peak pressure measured. The graph 
of Fig. B-10 plots this error for a family of wave durations for a linearly 
decaying pressure as a function of the time required for the shock wave to 
cross the gage's sensing area. 

CONCLUSIONS 

The conclusions to be drawn are based on the following considerations: 

1. The sensing element is circular with a diameter between 1/8 
and 1/4 in. 

2. The duration of any component wave is greater than 50 fisec 
for a triangular decay and 100 pisec for an exponential decay 

These conditions are not particularly restrictive and are met by the 
average modeled case encountered in the BRL experimental data. Based on these 
considerations, the following conclusions are made about the mapping of an 
incident pressure profile given as p(t): 

• The mapped pressure profile can be described by the relation 
p(t - R/U) and acts through the centerline of sensing area 

• If the mapped pressure profile is measured from the leading 
edge of the gage sensing surface, then it is in error R/U in 
time in relation to a fixed point on the model 

• The peak pressure can best be estimated by extrapolating the 
mapped pressure—time curve back in time R/U 

• The extrapolation of the curve provides a better estimate than 
use of the initial slope as the extrapolation model 
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MEASUREMENT 



CROSSING TIME, | (^sec) 


Fig. B-10. Error in Peak Pressure Reading Due to 
Gage Crossing-Time 




UR5 ■ 755-3 


• The error introduced by less than a 10% decay in the peak 
pressure across the gage element is not significant 

• The decay in the pressure profile created by the incidence 
of two rarefaction wave systems from different directions 
does not alter the previous conclusions. 
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Exhibit I 

FORCE VOLUME INTEGRATION OF A GENERAL WAVEFORM 


For the integration, the general waveform 


- r 

/ 

; x\ 

] r ( 

X \"1 

p L 1 

- c ( 

t - u) 

exp -ac ( 

‘-if)] 


( 1 . 1 ) 


has been selected. The terms are defined as: 

"p peak overpressure 

c inverse of the duration 
t time 

x distance 

U shock front velocity 

a Friedlander-type constant 

All definitions and nomenclature are consistent with the system described 
in the body of Appendix B. 


Substituting Eq. (1.1) into the force volume integral defined in the 
body of Appendix B gives 


2 x (y) 


v f = 2p exp (-act) J J £(1 - ct) exp + U X exp (FP)] 


dx dy (1.2) 


5 x x (y) 


If we define 


2p exp (-act) 9 


ac , c 

u ’ k - u 


0 = 1 - Ct 


B-22 


UR5I 755-3 


then Eq. (1.2) can be written as 

R x (y) 


= J7 

O x i (y) 


0 exp (wx) + kx exp (wx) I dx dy 


(1.3) 


Integrating with respect to (x) and rearranging gives 
R 


- /[($-*) 


X (y) 

2 


exp (wx) + — x exp (wx) 


dy 


x x (y) 


(1.4) 


0 k . , k . .. 

and if --—*■ = A and — = A , then 

w w 2 1 w 2 


R 


■ Jh 


x o(y) 


exp (wx) + A^x exp (wx) 


dy 


x (y) 
1 


(1.5) 


Evaluating the bracketed terms, 


A exp (wx)J 2 = A i exp |^ ex P ( w ^ R 2 -y 2 ) - ex p(-w^R 2 -y 2 )j 

X i 


= 2A 1 exp (wR) sinh w^ R 2 ~y 2 


( 1 . 6 ) 


and 


A g x exp (wx) 


2 = A exp (wR) R £exp ^w^ R 2 -y 2 ^ - exp ^-w ^ R 2 -y 2 ^J 


x 1 (y) 


•^R 2 -y 2 £exp ^w^R 2 -y 2 ^ + exp ^-w^ R 2 -y 2 ^ 


= 2A a exp (wR) ^R sinh w^ R 2 -y 2 + ^R 2 -y 2 cosh w ^ R 2 -y 2 ^ 


(1.7) 
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Equations (1.6) and (1.7) can be substituted into (1.5) to provide 


z = 2 exp 


R 

(wR) J* ^A i + A^R^ sinh ^ R 2 -y 2 ^ R 2 -y 2 coshdy 


( 1 . 8 ) 


Hyperbolic functions can be expanded in series (Ref. 3) as 


sinh x = 


x 3 x 5 x 7 

X + 3i + Ji + 7i + 


and 


cosh x = 


, X 2 X 4 X 6 

1 + 27 + 4i + sT + 


Combining these expansions with Eq. (1.8) and consolidating terms and letting 
A. + A R = A 

1 2 3 


z = 2 exp (wR) 


/ [(v + a 2 ) ( r2 ■ y2 ) 1/2 + 


1/2 /A g w 3 + 3 A 2 W 2 ’ 


3! 


(r 2 - y 2 ) 


3/2 


/A w 5 + 5A w 4 

U--— 

5! 


/ V 5/2 

(r 2 - y 2 ) + 


'A w 7 + 7A w 6 i 7/2 

3 2 1 / _ 0 o \ /* 


7! 


(R 2 - y 2 ) 


+ . . . 


dy 

(1.9) 


Integrating Eq. (1.9), evaluating at the limits, and rearranging terms provides 


p(t) = p exp 


H* -!)] 


/A w 3 + 3w 2 A 

( V + V + 2t2* - £ 


) R2+ ( 


A w 3 + 5 w 4 A > 

3 2 


4 ! 2 3 


R 4 


/Aw 7 + 7w e A \ 

+ ( 3 6!2 4 ) R6 + ••• 


(1.10) 
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which can be shown to equal 

T /, r M (, (. R\ . r 9 |"w[0 + k(wR + 3)] - k 

p(t) = p exp J^-ac^t ~ u)J J 1 " ^ ~ U/ + R |_ 8 

R 2 [w 3 (0 + wR + 5)] - k 
+ 192 

Approximate maximum values of the constants would be 
c = 10 4 sec” 1 

U = 10 3 feet per second 

R = 1CT 2 feet 

and assuming a value of a equal to 1, then 

k = w = 1 

For these values, the terms in the series of Eq, (1.11) after the first term 
are seen to add a negligible amount, and Eq. (I.11) can be written as 

p(t) = p exp £-ac(t - •§)][! - c(t - (1.12) 
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Exhibit II 
CENTROID LOCATION 


The centroid of the force volume for a triangularly decaying wave 
can be defined by use of a triple integral (Ref. 2). The centroid is needed 
to be known in the x,y plane. Since the y = 0 plane is a plane of symmetry, 
it contains the centroid (Ref. 5), thus the problem is to find the x coordinate 
x. The centroid x is defined by 



(II.1) 


The parameters used are consistent with those outlined in Exhibit I. 
integration yields 


x 



+ kR) 


+ 


kR] 

6 J 


V 


f 


The 


(II.2) 


where 0 = 1 - ct , and k = c/U. 


The force volume is described by 


7rR 2 


[p, + i <p 2 - p,>] = ^ <p, + p.) 


(II.3) 


and 


2pR 


x = 


(* + l kE ) 


P + P 
1 2 


(II.4) 


whe re 


p = p(l - ct) and p^ = p 


[> - • (‘ - ¥)] 
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Algebraically this reduces to 


R [(i ~ ct) + lf] 

[(1 - ct) + 


(II.5) 
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Exhibit III 

GENERAL CASE OF CENTROID OF FORCE VOLUME 

The centroid of the force volume (Figs. B-4 and B-6) in the x,y plane 
will determine the line of action in the z direction of the average pressure 
mapped by the transducer. If a plane of symmetry exists, then the centroid 
must lie on it (Ref. 5). For the cases being considered, the x,y plane for 
y = 0 is a plane of symmetry; thus one coordinate of the centroid is x,y = 0 
The x coordinate of the centroid x can be given by 


x V + RV 
w w c 



(III.l) 


as defined in Fig. B-8. Substituting in the values for V , V , and V in 

w’ c* f 

Eq. (III.l) gives 



(III.2) 


The worst possible case has x = 2R which reduces Eq. (III.2) to 

w 


x = R 


P 1 + (P 2 " P 1> 


• (^) 


(III.3) 


or 



(III.4) 
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Exhibit IV 

DECAY OF STRENGTH ACROSS THE SENSING ELEMENT 


Assume the peak pressure p varies with distance x across an element's 

K 

width of radius R according to the relation 


p k = p 


( l -s) 


(IV.1) 


where p is the peak incident pressure. Let a be the fractional spatial decay 
observed, and further assume no change in wave duration. A triangularly 
decaying wave would then be given by 


P(t) = P ( X - §)[* - c (* - §)] 


(IV.2) 


If we use the force-volume integral described in the discussion of the body 
of this report, then 


R x (y) 
>2 


R x (y) - 

2 cep 


V f - 2 


/ / P f 1 - ° (* - u)] dX dy - 2 / / 2R 

o X (y) o X (y) 


b - c (* • «)] 


dx dy 
(IV.3) 


By the techniques applied in Exhibit I the expresion for p(t), the average 
sensed pressure, can be shown to be 


p(t) = p { 1 - c 


M 1 -!) I]} 


(IV.4) 
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Appendix C 

PRESENTATION OF PRELIMINARY URS SHOCK TUNNEL TEST DATA 

The data on the air-blast-induced flows in full-scale rooms to be dis¬ 
cussed in this appendix were collected from tests in the URS shock tunnel. 

The aim of these tests was to evaluate the potential of the shock tunnel for 
aerodynamic testing and to recommend changes for more closely meeting future 
OCD test requirements. The findings reported herein are oriented toward the 
evaluative task of this study. An evaluation of the tunnel and of the sub¬ 
sequent improvements to the test procedures is being considered in another work 
unit and will not be included in this discussion. 

The room geometry of the two-dimensional test (the only one to be consid¬ 
ered herein) is shown in Fig.c-1. The geometry and wave form monitored at 
Station 1 were reproduced in the FLIC code as boundary and initial conditions. 

A direct comparison between full-scale experiment and numerical simulation 
could then be made. 

The first section of data to be presented is the pressure—time histories 
monitored in the room, the second section will consist of flow data. Varia¬ 
tions between the experimentally and the numerically derived data will be 
pointed out and then will be explained individually. It should be noted that 
the tests were run for tunnel evaluation purposes and, hence, the data are not 
representative of the results gained in present test series. Since limited 
full-scale data exist, and those data are not well documented generally, it 
was decided to make use of the available preliminary data from the tunnel. 

Data from two tests were used. The first test was the room shown in 
Fig. C-l, and the second was the same situation without a rear wall. The 
data presented for the downstream side of the front wall were taken from 
tests with no rear wall. The data for the front wall cannot be directly com¬ 
pared with the numerical simulations for the full time range recorded. The 
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Fig. C-l. Layout of Shock Tunnel Test 
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comparison becomes invalid at that point in time when the rear of the front 
wall senses the presence of the rear wall (approximately 2 traverse times of 
a sound wave across the chamber). 

In the second test, the input wave overpressure was lower, although a 
similar charge of primacord* was used in the compression chamber. The pres¬ 
sure-time records presented were adjusted accordingly by linear extrapolation. 
The normalizing factor used was the ratio of shock front overpressures monitored 
at Station 6. There is some error involved in using this simple adjustment to 
the same input condition; but, considering the accuracy of the data, it did not 
seem cost-effective to be any more elaborate. 

A moderately good agreement between numerical and experimentally produced 
pressure—histories can be noted at the input Stations 6, 7, and 11 (Figs. C-2 
and C-3). Variations later in the histories can be noted which correspond 
to a 10 to 20 percent lower experimentally determined pressure at fixed times. 
Excellent agreement was reached for the first 20 msec of each trace, though 
the values monitored were also lower in this range (generally about 10 percent). 
The absolute experimental accuracy was 10 to 15 percent; hence it is difficult 
to fully ascribe the cause of the variations to the code. 

Station 9, just inside the doorway, also showed good agreement between 
experiment and numerical simulation for about 20 msec. At about 30 msec the 
increase in pressure monitored in the experimental trace tapered off to a 
3 psi plateau. This measurement is not considered valid since an external 
overpressure in excess of 6.5 psi was monitored up to about 60 msec. Filling 


* 

Primacord is manufactured for wave-propagation characteristics; hence the 
primary manufacturing controls are not for equal charge weights per unit 
length. The resulting pressures developed by two different lots vary. 
Since the tests were 6 months apart, the primacord differed; hence equal 
quantities did not produce driven shock waves of equal strength. 
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Fig. 02. Pressure—Time Histories, Stations 6 and 11 
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STATION 7 



STATION 9 
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Fig. C-3. Pressure—Time Histories, Stations 7 and 9 
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of the chamber should have continued and the internal pressure should have 
increased to the external value. Otherwise a 3.5-psi discontinuity would 
have to exist between Stations 7 and 9 (a distance of 3.4 feet) and this 
is a physical impossibility. 

A similar problem was noted for Station Z on the rear wall. The moni¬ 
tored pressure stayed at about 4 psi after 35 msec. A break in the monitoring 
history can be noted at about 35 msec, which indicates a recording or monitor¬ 
ing system malfunction may have occurred at that time. 

Station N showed a much better agreement over the full range than did 
Station Z (Fig. C-4). The rarefaction trough in the 40-msec range was not 
reproduced as faithfully, and the characteristically lower experimental value 
again occurred. Also, the disparities in the pressure values yielded by the 
two methods are in the range of experimental error. If the input pressure 
history was lower (after 30 msec) in the experimental case, as indicated, the 
internal pressures must be accordingly reduced. 

The break in the monitored history, similar to that seen at Station N on 
the rear wall, was also seen in Stations H, I, K, and L on the rear of the front 
wall (Figs. C-5 and C-6). Stations H and K do show good agreement for about 
30 msec, however. The disruption of recording was traced to cable-noise prob¬ 
lems although the conditions previously noted also contributed. 

The two stagnation probes recorded for about 20 msec before being 
destroyed as illustrated in their pressure“time histories presented in 
Fig. 07. The stagnation values recorded (essentially the sum of dynamic and 
"side-on" pressures) compared well with the code for this period, an indica¬ 
tion that the code was predicting the flow as well as the pressure within at 
least a 15 percent value (assuming the experiment was correct). 

One major variation between the simulation and the experiment was the 
loss of the door jamb in the experiment. The high-intensity flows tore the 
wooden jamb from the steel frame, which caused a partial disruption of the 
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Fig. C-5. Pressure-Time Histories, Stations H and I 
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Fig. C-6. Pressure—Time Histories, Stations L and K 
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flow and may have been responsible for some of the deviations in the presented 
histories. The exact value of the error was not estimated because of the 
complexity of such a procedure and because, in the time frame used in the 
comparison, the perturbations produced by the error were relatively small. 

Flow visualization in the two-dimensional room was accomplished by means 
of ping pong balls and cloth flags placed throughout the room. A three- 
dimensional room in which a window replaced the door was also tested, but 
the data from it will not be discussed here. The data from the three- 
dimensional room were superior to those from the two-dimensional case and 
were used in establishing the data reduction procedures. They were also used 
in the data reduction process to check for errors, because similar flows 
occurred in both situations. 

The raw data were in the form of a series of pictures taken with two 
cameras. The location of the cameras and the ping pong balls are shown in 
Fig. C-8. The arrival of the shock wave was monitored by the movement of 
cloth flags along the roof of the room. Distances were measured by means of 
vertical and horizontal grids on each wall. In the film plane, the balls 
appear against the background of the grids. The true spatial positions of 
the balls were established from the geometry of the two camera projections 
(congruent triangles). 

The tests were originally planned as preliminary trials for another OCD 
work unit, as previously discussed. The data of use to the computer evalua¬ 
tion effort of this work unit were removed and adapted for present use. The 
data reduction and evaluation techniques being developed for that work unit 
were used in their preliminary form with some adaptations. 

The first test was the two-dimensional configuration. Lighting 
difficulties were encountered and were corrected for the second test, the 
three-dimensional window configuration. Since the numerical simulation eval¬ 
uations needed two-dimensional data, the poorer test had to be used. The 
motion of only three balls could be continuously followed for a sufficient 
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period to provide the raw data to establish their trajectories. Additional 
qualitative data were, however, evolved for short time spans during the flow 
history. Though useful in developing the conceptual model of the flow, the 
short-time-span data will not be used to develop quantitative information on 
the flow. 

The trajectories of ping pong balls, though indicating the flow direction 
and magnitude, are not a direct measure of the flow velocity. As a guide to 
understanding the underlying process, the question of how rapidly a ping pong 
ball can be accelerated by the n jet” flows was considered. The acceleration 
determines what fraction of the particle velocity of the flow the ball will 
attain as a function of time. These acceleration curves will provide some 
insight into what the magnitude of the driving flows must be in order to 
obtain the monitored ball translation velocities. 

The physical model of the acceleration process is simply a Newtonian 
force balance on the target ball. The force is provided by the drag of the 
flow field surrounding the ball, i.e., the flow velocity experienced by the 
ball is relative to itself, not to laboratory coordinates. The entire system 
can be transformed to the fixed laboratory coordinates simply by defining the 
driving flow in terms of the velocity differential between the flow and the 
sphere. The system is presented in Fig. C-9. 

Equation (C.3) is nonlinear both in x and t (see Fig. C-9). For 
illustrative purpose, assume a step wave is incident on a target ball. 
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p(x, t) 
p(x,t) 
u(x,t) 






x coordinate 


p(x, t) 
p(x, t) 
u(x, t) 






x coordinate 


F 


p(x, t) 


[u(x, t) 



(C.l) 


F 



(C. 2) 


dx 

dt 


r = \ p(x,t) [u(x, t) - x] ! 


V 

m 


(C. 3) 


Fig. C-9. A Schematic of the Sphere Flow Field Interaction 


For this case, Eq. (C.3) can be written (for all t ^ 0) as 

/c a\ 


d(x) 


= \ p C u " k T 


D 

m 


(C. 4 ) 


assuming is a constant. Ahlers (Ref. 36) considered Eq. (C. 4) to be of 

the Ricatti type and solved it by a transformation of coordinates. For the 
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step conditions assumed, Eq. (C.4) can be directly integrated without any subs 

titutions. Using either method, the solution to Eq. (C.4) can be shown to be 
(using Ahler's nomenclature) 


x = u + 


C 3 /3 


1 + C 3 t 


and 


x = C 4 + ut + )3 i (1 + C,t) 


where C. and C 4 are constants of integration and 


P = 


-2 


AC. 


D 


m 


The value of r equals -1 if x ^ u and equals +1 if x ^ u. 


(C. 5 ) 


(C. 6 ) 


(C. 7 ) 


In his work, Ahlers further considered the case in which p u, and p 
were not invariant with x and t. The previous solution is invalid for 
spatial and temporal variations in thermodynamic and flow properties. 

Ahlers' approach was to represent the above solution in a difference format. 
The differencing steps in this format could be selected so that average values 
of p and u over the time step could be used without introducing significant 
errors. The following equations present Ahlers' difference equations with 
errors in the original work corrected. 


x 

m + 


1 


(* 


U + £ 
m p 


m 


u ) (3 
m m 


(x - u ) At 
m m 


m 


(C.8) 


x 

m + 


1 


x + u At + j3 i 
mm m n 



(C. 9) 


One added variable, which is part of the term 0 , is the drag coefficient 
The values of the drag coefficient for a sphere are presented in Fig. C-10 
as a function of flow Mach number. As the space-time variations in the flow 
occur, the corresponding corrections in the drag coefficient must be made. 
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Fig. C-10. Drag Coefficient for a Sphere 
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A branch in the curve can be noted at a Mach number of approximately 0.775. 

The branching is dependent upon whether turbulent or laminar flow occurs. The 

5 

transition Reynolds number is 4 x 10 . If the inflow velocities range from 

100 to 1000 fps at standard sea level densities and temperatures, then R 

„ n 

4 5 

would vary from 7.8 x 10 to 7.8 x 10 , which includes the transition between 
laminar and turbulent flow. 

Before making the comparison of numerical and experimental simulations, 
the original consideration of the acceleration time of the ping pong balls 
to a step drag profile will be discussed. The profile selected will be for 
a jet of fluid with the conditions described in Ref. 2 for free-flow jets 
initiated by a shock wave impinging on a chamber. From these values, the 
time it takes a ping pong ball to reach a percentage of the flow speed (U) can 
be calculated using Eq. (C.9) and the drag coefficients presented in Fig. C-10. 
A constant density and viscosity were assumed corresponding to a fully 
expanded flow. The drag coefficients were based on calculated values of 
Reynolds and Mach numbers for the flow in each time step. The computations 
were made using a short computer program which selected the proper drag 
coefficient based on Mach and Reynolds numbers and then used these values to 
calculate velocity data. The results are presented in the curves in Fig. OIL 

These curves provide some interesting insights into the motion of the 
balls. The rather long acceleration times and the asymptotic behavior of the 
balls' velocity—time profile reflects the decreasing driving force as the 
balls approach flow velocity. At 10 msec, the balls are at 25 to 40 percent 
of the free-stream particle velocity; at 100 msec, the balls are 75 to 85 
percent of the particle velocity; and at 1 sec about 97 percent of the particle 
velocity. The crossing of the curves for AP = 9 and 12 psi and some of the 
inflections in the curves correspond to changes in the drag coefficient 
created by transitions from laminar to turbulent flow as defined by the flow 
Mach number and Reynolds number. Though not a precise prediction procedure, 
the test for the drag conditions is considered sufficiently accurate for the 
exploratory and explanatory purposes of this study. 

In reviewing the equations and the curves for the drag coefficient, the 
question arises: If the trajectory of the target object is known, can the 
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flow field be inferred from it? The drag coefficient, as presented in 

Fig- C-10, is a function of Mach number and Reynolds number. The Mach and 

Reynolds numbers are, in turn, functions of the velocity differences between 

the flow and the object, the viscosity and density of the gas, and the object's 

diameter. The experimental data provides x x , and At from the ball's 

m+l m 

trajectory. 

More explicitly, an expression for u , the flow velocity during the 
period At, can be derived from Eq. (C.8) in the form 

At 


m+l 


X - (x - u ) (x 
m m+l m m 


u ) 
m (3 


0 


(C.10) 


m 


which reduces to 


- C, u 
m 1 m 


c sA 


(C.ll) 


where C 1 , C 2 , and C 3 are constants determined from the initial conditions, 8 

m 

is determined from these conditions and C^, the drag coefficient, which can be 
expressed in functional form as 


D 


= f (u 


m 


a > P. u) 

m m m 


(C.12 ) 


If the flow conditions at one point in the flow field are known and a 
single condition is known at another point of interest then Expression (C.10) 
could be solved by iteration. To visualize this, consider a one-dimensional, 
constant-area, homentropic flow field. For these conditions, the isentropic 
relationship and the conservation of energy and mass are given by 
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p u = constant = /2 D 0 (C.15) 
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These equations can be combined to yield the expression 

(C.16) 

which cannot be solved explicitly. 

Assume p is known. Then the conditions needed to determine sound speed 
m 1 

and viscosity can be derived from Eq. (C.13) and the equation of state. Under 
these conditions Eq. (C.ll) could be rewritten as 

u 2 - C u + C f (u ) = CL (C.17) 

m l m 4 m J 

where C is a constant and C = f(u ). Even with these "most favorable condi- 
4 Dm 

tions", an iterative solution is required since u must be known to determine 

m 

f(u ). The conclusion to be reached is that the procedure is somewhat less 
m 

than an optimum method of directly developing flow speed measurements. In 
reality, it is an approximation if only experimental data are used. Further 
complicating the process is the fact that arguments can be posed as to the 
true values of transient drag coefficients and the effects of lift and spin 
on the trajectories of the balls. 

In the present effort some of the above problems were avoided by putting 
the target objects into the flow field simulated by the FLIC code. The flow 
properties are precisely known in the computer simulation; hence if we assume 
steady-state drag coefficients apply and that no lift or spin effects occur, 
the trajectories of objects in the simulation can be determined. The 
calculated trajectories can then be compared to the experimental trajectories. 

In the calculation procedure the thermodynamic and aerodynamic parameters 
of the flow field were established via interpolations for the space-time 
points in the flow. The calculation procedure was simplified by taking very 
small difference steps and using constant values across it for all parameters. 
The simplified Newtonian equations were placed in a computer program that 
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also included subroutines for calculating the drag coefficient, viscosity, 
Reynolds number, sound speed, and Mach number each as a function of the 
input variables. The program was checked against the more complex one 
used to calculate the acceleration data. The agreement between the data 
evolved from both programs was found to be consistent within a 1% error. 

The results of the computations are presented in Figs. C-12 through 
C-14. Fig- C-12 presents the trajectory data. The variation between 
numerical and experimentally determined trajectories ranges from 1 to 3 
msec for fixed points on the trajectory curves. This variation corresponds 
to one-half to two frames of the film, and also corresponds to about the 
best zero-time measurement (shock arrival) for the method used. The agree¬ 
ment was surprisingly good. If, for example, a 1- to 2-msec error in shock 
arrival time did occur, as is indicated by the velocity curves of Figs. C-13 
and C-14, then on the adjusted time scale the numerical and experimental 
trajectories superimpose. 

The errors in the translational velocity histories shown in Figs. C-13 
and C-14 were larger. The velocities of the experimental objects were 
initially low. They were then subject to some unrealistically large accelera¬ 
tions and became slightly greater than the numerically calculated values. The 
variation in the slope between experimental points caused these errors, as 
will be shown in succeeding paragraphs. Before discussing the errors, the 
effect of smoothing the data will be reviewed. 

The smoothing technique used was a five-point Lagrangian interpolation 
which was based on methods presented in Ref. 38. The interpolation routine 
was a subroutine of the program being developed to reduce data from object 
translation for Work Unit 1154G. Of the available data, only ball No. 4 had 
sufficient points to warrant applying the smoothing techniques. The results 
are presented in Fig. C-15. 

The experimental curve drawn from the smoothed data has a noticeably 
different profile than the curve for the unsmoothed data. The inflections 
correspond to some of the rarefaction wave motions, but appear to have an 
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Fig. C-13. Comparison of Velocity Time Histories of a Ball in 
Shock Tunnel Tests to Numerical Simulations 
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Fig. C-14. Comparison of Velocity-Time History of Balls in 
Shock Tunnel Tests to Numerical Simulation 
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overly strong effect. The numerical simulations do not show the same 
pronounced decrease in acceleration, which would indicate that the simulation 
spread the steep rarefaction behind the shock wave due to the limited number 
of cells in the entry area. Although this may occur, the smoothed data 
appeared to have an inconsistency sufficiently serious to question their 
validity. The smoothed curve shows the ball decelerating, indicating a flow 
under 45 fps, well below those observed in either the numerical simulations 
or the BRL experiments (Refs. 4 and 5). 

The accelerations in some sections of the smoothed curve appear to be 
too great in view of the acceleration curves presented earlier. For example, 
at 110 fps, the smoothed curve showed the ball accelerated to 210 fps in 2 
msec.. This acceleration can only be approached in one small section of the 
AP = 9 psi curve in Fig. C-ll. This result would presume an instantaneous 
change in flow properties at 20 msec. These arguments indicate that the 
smoothing improved the presentation, but did not eliminate all the experimental 
variation. The conclusions drawn are that the smoothed data leveraged actual 
changes in the flow and that a more extensive analysis is warranted in future 
work to determine what smoothing techniques should be used. 

Even with this variation, the final velocities were within 10 percent of 
each other. Part of this is attributable to the effect of experimental error. 

A review of the equations will yield the presence of a second effect which 
will partly explain why deviations would tend not to be maintained. The 
driving force is directly proportional to the square of the object-to-flow 
field velocity differential. A slower object speed relative to a fixed flow 
speed will result in the object being accelerated faster. The more rapid 
acceleration will reduce the velocity difference between objects and the 
flow stream. If the same object is subjected to identical flows in two 
simulations but has a slower initial velocity in one simulation, then for 
the above reason the slower object will accelerate more rapidly reducing the 
object velocity differential between simulations. 
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In a sense errors are somewhat self-limiting, unless the spatial- 
temporal relationships between simulations have a large variance. The 
trajectories result from different driving forces which, in turn, increase 
the error, as was the situation for smoke particles in the BRL experiments. 

The velocities of the balls did not become very large (about 200 fps in the 
time frame of interest). The spatial variation for the smaller displacements 
and the accompanying errors were much smaller. The shifts noted in the scale 
model comparisons would have occurred on a reduced scale. In that the spatial 
variation between experimental and numerical simulations did not occur, i.e., 
the object trajectories compared reasonably well, it is concluded that the 
variations in velocities were basically attributable to experimental errors. 

On this basis the code-to-experiment dynamic pressure correlations would be 
judged to be within 10 percent. 

To develop a better understanding of the experimental data and the 
meaning of the experimental errors a brief error analysis will be undertaken. 
Fig. C-16 presents the nomenclature for the x,t data to be considered. Each 
x,t set corresponds to a data point reduced from the films taken during the 
experiments. 

The data sets correspond to the discrete points in monotonically 

increasing time from which the positibn-time curves presented earlier were 

drawn. The mean value theorem states (Ref. 30), "Let f be continuous on 

fa,b] and let f'(t) exist for a < t < b. Then, a point t can be found, with 

o 

a < t < b, such that f(b) = f(a) + (b-a)f f (t )." If the position time curve 
o o 

is defined as f, and if the displacement in a time interval At (given by b-a) 

is f(b) - f(a) then a velocity exists at some time in the interval At which 

corresponds to [f(b)-f(a)]/(b-a)." The velocity so defined is simply the change 

in displacement (Ax) divided by the time increment (At) . This argument states 

that it occurs at some time t in the interval At, but this does not mean that 

o 

it must be the average velocity in that interval, nor does it locate the time 

t . 
o 
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In taking this "average" value and ascribing it to the mid-point of the 
interval At an error results. An exact estimate of the error involved in this 
procedure is dependent upon the function f in the interval At. Crandall 
(Ref. 40) provides several methods of analyzing the error for parabolic 
curves through three adjacent points. Similar arguments are presented by 
Hamming (Ref. 41). Both presentations show the error to be a function of 
smoothness of the curve and the step size, i.e., the closer the segment of 
the curve is to being linear the smaller the error and the nearer the average 
value falls to the mid point of the interval At. 


The trajectory curves shown in Fig. C-12 can be quite accurately 
represented by straight line segments for the intervals shown. The indication 
is that this is not the major source of error contributing to the observed 
deviation in the velocity-time profiles. The source of error can, however, be 
shown to be related to small measurement errors in the x-t sets of the 
trajectory data. 


By using either the law of the mean or the definition of a derivative 
we can define* the velocity under the above conditions as 


dx _ Ax _ f(a) - f(b) 
dt “ At a - b 


(C. 18 ) 


which can be rewritten for the nomenclature in Fig. C-16 as 



(C.19) 


For convenience sake the set of points x ,t is considered to be between the 

n n 

two experimental data sets n + 1 and n - 1. If we assume displacement 


*This is by definition only, for in reality 


dx 

dt 


lim 
Ax — 0 
At—0 


Ax 

At 


lim 

(a-b) 


f(a) - f(b) . 
a-b 5 


If x 


f (t) 
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measurement errors (e ) occur, then the true value of (Ax/At) can be given 

x n 

as 



(x 


n+1 


- x 


n-1' 


n+1 


t 


n-1 


(C.20) 


An experimental measurement error of ±0.5 inches will be assumed 
which is not completely unrealistic. The data had a At value of 1.7 msec 
and corresponding Ax’s in the range of 3 to 6 inches (At and Ax correspond to 
n). The error in the computed velocity would be given by the relationship 


/Ax\ _ _.+ fx 

U*7 n VAt/ n _ :"At 


(C.21) 


The percent error 


(e ) would be given by 
vp 



( 100 ) 



(C.22) 


(C.23 ) 


1 

± + 1 


( 100 ) 


(C.24) 


For the conditions Ax = 5 in. and e = 0.5 in., the error would be ±9.9 

x 

percent. The geometry involved in the reduction process of the projections 
includes measurements up to 20 ft; therefore, a 0.5 in. error corresponds to 
a 0.2 percent error. The ping pong balls also move during each frame and 
estimates of their centers must be made from blurred images. Errors of 1/8 to 
1/2 in. are not therefore unreasonable. A more thorough analysis of these 
errors is being presently undertaken in Work Unit 1154G. 
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A second source of error could be the variations in the framing speed. 
These would be slower variations and most likely would follow a cyclic or 
smooth curve. The smoothed data in Fig. 015 showed this form of variation. 
An approximate sinusoid with about a 7-msec natural period (150 Hz) appeared 
This frequency is approximately the frequency of a harmonic which appears in 
the electrical power at the experimental facility. 


The true value including errors in both timing and displacement measure¬ 
ments would be given by the following relationships. 


'Ax\ 


Ax ± € 
_? 

At ± € , 


Ax 


At ± 


At ± €. 


(C.25) 


The percent error in the measurement would be given by 


€ 


vpt 


(Ax/At±e t ) ± (e^/At' x e t ) - (Ax/At) 
(Ax/At+e t > ± (e x /At±e t ) 


( 100 ) 


(C. 26 ) 


1 * <£ t /e I ) V At) n 

1 ± (Ax/e ) 

x 


( 100 ) 


(C.27) 


If e t “*~0, then Eq. (C.27) reduces to Eq. (C.24). For the same conditions 
as previously discussed (Ax/G = 10), and with 6 /At = +0.1, the maximum total 

X V 

error would be 22.2 percent. For e /At = +0.02, the maximum total error would 
be 13.3 percent. The timing error could therefore significantly contribute 
to the total measurement error. 


The preceding analysis shows how small measurement errors in trajectory 
data can be compounded by the conversion of the data to velocities. For this 
reason, the spatial trajectory data was used in the comparison of numerical 
to experimental simulations. The disadvantages and problems of computing 
translations, discussed earlier in this appendix, would cause smaller errors 
than differentiating the trajectory curve. Using the trajectory comparison, 
an estimate of the correlation between numerical and experimental dynamic 
pressure data would appear to be better than 90 percent. A more precise 
estimate of variation would require a more thorough error analysis of the 
experimental data. 


C-31 


Copies 


DISTRIBUTION LIST 


45 Office of Civil Defense, Office of the Secretary of the Army, Attn: Research Administration Office, Pentagon, Washington, D.C. 20310 

1 Army Library, 1A 518, Pentagon, Washington, D.C, 20310 

1 Assistant Secretary of the Army (R&D) Attn: Assistant for Research, Washington, D.C. 20310 

1 Chief of Naval Research, Department of the Navy, Washington, D.C, 20360 

1 Commander, Naval Supply Systems Command (Code 0611C1), Department of the Navy, Washington, D.C. 20360 

1 Commander, Naval Facilities Engineering Command, Research and Development (Code 0322C), Department of the Navy, Washington, D.C. 20390 

1 Mr. Richard Park, Advisory Committee on Civil Defense, National Academy of Sciences, 2101 Constitution Avenue, N.W., Washington, D.C. 20418 
20 Defense Documentation Center, Cameron Station, Alexandria, Virginia 22314 
1 Mrs. Joanne S. Gailar, Civil Defense Research Project, Oak Ridge National Laboratory, P. 0. Box X, Oak Ridge, Tennessee 37830 

1 Mr. Norward A. Meador, Shelter Research Division, Office of Civil Defense, Dept, of the Army - OSA, Washington, D.C. 20310 

1 Chief of Naval Personnel, (Code Pers M 12), Department of the Navy, Washington, D.C. 20360 
1 u. S. Naval Civil Engineering Laboratory, Port Hueneme, California 93041 

1 Director, Disaster and Defense Services Staff, Agricultural Stabilization and Conservation Service, U.S. Dept, of Agriculture, Washington, D.C. 20250 

1 Mr. Bill Miller, Department of Civil Engineering, 307 More Hall, University of Washington, Seattle, Washington 98105 

1 Mr. Carl Koontz, Department of Civil Engineering, Worcester Polytechnic Institute, Worcester, Massachusetts, 01609 

1 Mr. Robert Bailey, School of Civil Engineering, Civil Engineering Building, Purdue University, Lafayette, Indiana 47907 
1 Mr. John A. Samuel, Department of Mechanical Engineering, University of Florida, Gainesville, Florida 32601 
1 Mr. G.K. Vetter, School of Architecture, University of Colorado, Boulder, Colorado 80302 

1 Miss Nancy K. Barberii, OCD Professional Advisory Service Center, University of Arizona, Tucson, Arizona 85721 

1 Mr. Dick Rummer, 101 Eng. A, Pennsylvania State University, University Park, Pennsylvania 16802 
1 Director, Defense Atomic Support Agency, Attn: Jack R. Kelso, Washington, D.C, 20301 
1 Director of Research and Development, Office of Emergency Preparedness, Washington, D.C. 20504 

1 Director, Civil Effects Branch, Division of Biology and Medicine, Atomic Energy Commission, Attn: Mr. L.J, Deal, Washington, D.C. 20545 

1 Mr. George N. Sisson, Director, Shelter Research Division, Office of Civil Defense, Dept, of the Army - OSA, Washington, D.C. 20310 

2 Air Force Special Weapons Laboratory, Attn: Technical Library, Kirtland Air Force Base, New Mexico 87117 
1 Los Alamos Scientific Laboratory, Attn: Document Library, Los Alamos, New Mexico 87544 

1 Chief of Engineers, Department of the Army, Attn: ENGME-RD, Washington, D.C.20315 

1 Chief, Joint Civil Defense Support Group, Office, Chief of Engineers, Department of the Army, Gravelly Point, Washington, D.C. 20315 

1 Director, Army Materials and Mechanics Research Center, Attn: Technical Library, Watertown, Massachusetts 02172 

1 Director, U.S. Army Ballistic Research Laboratory, Attn: Document Library, Aberdeen Proving Ground, Maryland 21005 

1 Director, U.S. Army Ballistic Research Laboratory, Attn: Mr. William Taylor, G. Coulter, Aberdeen Proving Ground, Maryland 21005 

1 Director, Defense Atomic Support Agency, Attn: Technical Library, Washington, D.C. 20301 

1 Director, U.S. Army Engineer Waterways Experiment Station P. 0. Box 631, Attn: Nuclear Weapons Effects Branch, Vicksburg, Mississippi, 39180 

1 Director, U.S. Army Engineer Waterways Experiment Station P.0. Box 633, Attn: Document Library, Vicksburg, Mississippi 39180 

1 District Engineer, U.S. Army Engineer District, Omaha, Attn: Chief, Engineering Division, 6012 U.S. Post Office and Courthouse, Omaha, Nebraska 68101 
1 Mr. Carl K. Wiehle, Civil Defense Technical Office, Stanford Research Institute, Menlo Park, California 94025 
5 Mr. William L. White, Civil Defense Technical Office, Stanford Research Institute, Menlo Park, California 94025 

1 Mr. Werner Weber, Director, New York State Civil Defense Commission, Public Security Bldg., State Office Bldg. Campus, Albany, New York,12226 
1 Agbabin-Jacobsen Associates, 8943 South Sepulveda Boulevard, Los Angeles, California 90045 
1 Amman and Whitney, 111 Eighth Avenue, New York, New York 10011 

1 Mr. Arthur D. Caster, Chairman, Coordinating Committee on Civil Defense, American Society of Engineers, 2864 McFarlan Park Dr.,Cincinnati, Ohio 45211 

1 The Dikewood Corporation, 1009 Bardbury Drive, S.E., University Research Park, Albuquerque, New Mexico 87106 

1 General American Transportation Corporation, General American Research Division, 7449 North Natchez Avenue, Niles, Illinois 60648 

1 Husdson Institute, Quaker Ridge Road, Croton-on-Hudson, New York 10520 

1 Bell Telephone Laboratories, Inc., Attn: Mr. R. W. Mayo, Whippany Road, Whippany, New Jersey 07981 
1 Dr, Eugene Sevin, IIT Research Institute, 10 West 35th Street, Chicago, Illinois 60616 

1 Dr. Harold Brode, The RAND Corporation, 1700 Main Street, Santa Monica, California 90401 

1 Research Triangle Institute, P.0. Box 12194, Research Triangle Park, North Carolina 27709 

1 Mr. Luke J. Vortman, Division 5412, Sandia Corporation, Box 5800, Sandia Base, Albuquerque, New Mexico 87115 
1 URS Research Company, 1811 Trousdale Drive, Burlingame, California 94011 

1 Massachusetts Institute of Technology, Department of Civil and Sanitary Engineering, Cambridge, Massachusetts 02138 
1 Dr. Nathan M. Newmark, University of Illinois, 111 Talbot Laboratory, Urbana, Illinois 61801 

1 Dr. William Hall, University of Illinois, 111 Talbot Laboratory, Urbana, Illinois 61801 

1 Dr. Merit P. White, University of Massachusetts, School of Engineering, Amherst, Massachusetts 01002 

1 Dr. Abner Sachs, Institute for Defense Analyses, 400 Army-Navy Drive, Arlington, Virginia 22202 

1 The Vertex Corporation, 10400 Connecticut Avenue, Kensington, Maryland 20795 

1 Chief of Engineers, Department of the Army, Attn: ENGMC-EM, Washington, D.C. 20315 

1 Dr. C. S. White, President-Director, Lovelace Foundation, 5200 Gibson Boulevard, S.E. Albuquerque, New Mexico 87108 

1 Dr. Charles Osterberg, Acting Chief, Environmental Sciences Branch, Division of Biology and Medicine, U.S. Atomic Energy Commission, Washington, D.C. 20545 

1 Mr. J. J. Davis, Effects Evaluation Division, Nevada Operations Office, U.S. Atomic Energy Commission, Las Vegas, Nevada 89101 

1 Mr. Eugene F. Witt, Bell Telephone Laboratories, Inc., Whippany Road, Whippany, New Jersey 07981 

1 Mr. Paul Zigman, Environmental Science Associates, 770 Airport Boulevard, Burlingame, California 94010 



Unclassified 
Security Classification 


DOCUMENT CONTROL DATA - R & D 

_ (Security classification of title, body oi abstract and indexing annotation must be entered when the overall report is classified) 


1. ORIGINATING ACTIVITY (Corporate author) 


URS Research Company 
1811 Trousdale Drive 
Burlingame. California 94010 


2«. REPORT SECURITY CLASSIFICATION 

_ Unclassified _ 

26. GROUP 


3. REPORT TITLE 


The Air-Blast-Induced Environment Within Civil Defense Blast-Slanted Shelters 


4. DESCRIPTIVE NOTES (Type of report and inclusive dates) 


5 AUTHOR(S) (First name, middle initial, last name) 


Joseph F. Melichar 


6. REPORT DATE 

November 1969 


8a. CONTRACT OR GRANT NO. 

DAHC20-67-C-0136 

6. P ROJEC T NO. 


c. OCD Work Unit 1123F 


10. DISTRIBUTION STATEMENT 


7a. TOTAL NO. OF PAGES 

236 

7b. NO. OF REFS 

68 

9 a. ORIGINATOR'S REPORT N UMT3 E R(S) 

URS 755-3 


I 9b. OTHER REPORT NO(S> (Any other numbers that may be assigned \ 

this report) 



This document has been approved for public release and sale; its distribution 
is unlimited. 



12. SPONSORING MILITARY ACTIVITY 


Office of Civil Defense 


13. ABSTRACT 

The study evaluated the potential usefulness of two dimensional hydrodynamic computer 
codes to the civil defense shelter research program, and also continued work on 
evolving a conceptual model of air blast-induced flows in "blast-slanted" basement 
shelters. The review of the numerical simulation techniques included the theoretical 
aspects of the code usage as well as questions of practical usefulness. The 
numerical simulations were compared against experimental simulations and it was 
concluded that they were less expensive and more flexible. It was also included 
that a blend of both simulation techniques is necessary to the shelter research 
program and that the blend between the two techniques and analysis is directly 
related to specific research goals. The numerical simulations, experimental simu¬ 
lations, and analytical techniques were used to extend the understanding of basic 
flow process of the air blast transmission process into the shelter. These results, 
coupled with past work, were used to show the economic infeasibility of developing 
a generalized conceptual model of the flow into blast-slanted basement shelters. 
Questions were raised as to the need for such a model and recommendations as to 
future research needs were made. 


DD 


OftM 

• MOV #S 


1473 


MFLACII DO FORM I47». 1 JAN S4, WHICH IS 
OtSOLKTK FOR ARMY USB. 


Unclassified 

Security Classification 

















Unclassifi 


Security Classification 



Unclassified 


Security Classification 











IS 3 o 
£ 


V w 
r ft w 
3 ft P ffl 

3d) cp 
J P X P P 

i c <d s 

< « T3 <H 13 

> P 3 Q) e 


I H ffl OP 

ft I ffl P p 

is I 3 P rt 

O' O P ■ 

C P BJ 3 I 

P C P 6 

ft -H 

wo w 


I P C W P ■ 

3 O C C 

* ft P o 0) 

QJPP E < 
, O rt P P ‘ 

3 P w P 


<D ft W P 

P W X P ffl 

P P ffl ffl P 

(D W 'OP 1 

ft >> •> C 0) 

w h to 3 ft • 

a C w 

0) C O 0) 

ft d -H £ ffl • 

p p p ft 

ft a p 

O C P -o 

P rt 3 C O 

e ffl p 


ffl £ O P ffl 

O ft P O 

ffl O P ft 0 

C ffl ffl ffl P , 

p 6 w ft 
w 3 3 

•HOC c 


? ft ffl I 
P ffl O P 
rt p O ffl 
3 > ■* 

p ffl ffl 

ft p ft wj 


■ C ffl OP 

« O ft ft ffl w 

p o p rt 

P *o P P 

rt C rt h 

H ffl ffl rt 

3 p w o p 

Eft ffl P p 

p P p rt 


ffl ffl O £ 
P *o P 
ffl o w 
? £ w -O 
rt ffl 

* rH <D 

X rt -3 C 
P 3 ffl 
O P w xt 
S ft P O 
ffl rt P 


ffl w 

z ft w 

13 -O P ffl 

ffl ffl C 

O 3 P 13 P • 

•rt £ £ ffl 3 

g P rt 13 P 1 


s 13 p rt 

COP 

C p rt 3 
P C p E 
,ft ft p 
wo w 
■ £ ffl p 

OP 0 p I 


•H o 13 P 0 ffl 

13 P ffl p p E f 

ft o rt p p ‘ 

0 3 p w P 

£ ft 13 3 ffl ffl 

P O C E 3 ft 


3 WO ffl ffl P 

J5 hfl H 

P ffl p p rt ft 

rt w ffl WE 

p C ft P 3 O 

p ffl o o o 

CPS ffl 

ffl ffl £ 13 ffl I 
P 13 P ffl O P 

O rt p o ffl 

ft p 3 > * 


ft rt o P 

■> ft ffl ft o 

I <H p ft c 


rt p o w o 

P 13 C rt g 
3 C O 0 P 


c o p ffl 

ftp O * s* P 

O P 13 0 0 O 0 

ffl ffl ffl P ft P P 

p e w ft w p 

3 3 ft • 

O C C O ffl ffl ffl 

£ ffl O p ft ffl ft 

p ffl p p p C rt 

ft ffl w -o s 

ffl H ? w ffl p ffl 

ft p w O ft ffl 

p w e 3 p p 


I (DP O’ rt P 13 P 

i ffl rt p p ffl o w 

’^OCP^BWft 

p ho ft rt ffl 

! ffl o P - P ffl 
iftftfflwftrtUC 
i o p rt P 3 ffl 

'•OP P O P W ft 
iCrtPftftftPO 
i ffl ffl rt ffl rt P 

ipw oppoPrt 

! ft ffl P P W C ffl 


> P rt ft ft 13 5 

P C P P ffl ffl 

> P rt P N W P 

: o p * p c 3 
i ffl ft 0 pop 

1 ft C 13 rt P 3 


rt • P 3 

WOE 
fcJD P ffl P , 

C ffl ft w 

P p W 


O ft p o • 
> w rt p 

ffl O P r 
P -rt ffl 
C C +» E 
O ffl ffl 3 
E P C 
I X ffl o 
P W ffl ffl 

o rt ft ft 

S ft p Eh 


ft P P p w 

ffl rt o rt rt 

13 ffl ffl p ft 

3 w p c 

p ffl P ffl P 
, o p ft E O 


ft C ffl o 
3 ffl 3 P 
0 bfi O’ 


P ffl ft 13 
> w C 
ffl ffl ffl 
w 13 P E 
ffl C E 
ft P ffl O 
Eh O E O 


ft ffl P P 

w n in ffl c 

I I ffl ft d 

X ft CS E 

P ft ffl J* 

<C O CO > P 

os o o 

W W ft 55 S 


rt rt • p 3 

3 0 WOE 
P p bfi P ffl P , 

rt C ffl ft w ■ 

> W P p w 

ffl ffl > p rt p ’ 


w ffl c c p g 

p O ffl ffl 3 

ffl 3 E p c| 

ft ft ft ffl O 

Eh E P W ffl ffl 

o o rt ft ft 

O » ft P Eh 


CO - ffl A fi ffl O 
rt P WO 3 ffl 3 P 
„ £ O O bC O’ 

S 13 O P O W 

rt ffl p ft rt rt 
p p p * • 

ho rt rt £ w bowlw 

0 P p OP C PC 

P ffl 3 p P P (D O 

ftp S «H 3 ft P P 


P P P W > (A C 

rt o rt rt ffl ffl ffl 

ffl ffl P ft W ft P E 

w P C ffl C 1 
(DP ffl P ft P ffl o 
PftSOHOBO 


ft (DPP 
CQ CO W ffl c 

1 I ffl ft ft 

OS 1C OJ E 

p in ffl ft 

< t* w > p 

os o o 

WWDS5J 
ffi OS 
Eh ft 


r ft w 

3 13 P ffl 

D ffl CPI 
3 p 13 P P C 

3 C ffl 3 P 

h rt ft p xi 

■> p 3 ffl C 0 


3 P W OP rt 

1 ft! ffl p p & 

Jr I 3 p rt 

o' o P p 
J C p rt 3 p 

: p C p E 

j ft ft p 

wo w * 

- s ffl p ffl 

sop o p p 

i p rt ft 

i P C W p P 

I O C C X 

> ft P O ffl ffl 

I ffl P p E p 

l o rt p P p 

3 p w P 

! ft 3 ffl ffl ffl 

> C E 3 ft P 

I P p O’ X O 

I I w ffl E 


w p « ft rt 

P P ffl p ■ p ft 

ffl ft W p w 

P w x p <D rt 13 

p p ffl ffl p ffl ffl 

ffl W x! P P P 

ft >> -, C ffl c c 

w P w 3 ft p rt 

rt c w p 

ffl £ O ffl o w 


ft rt p o w 

o C p ft C rt 

P rt 3 C O O p 

E ffl P O ft| 

>. w p p c ffl 

P ffl w X p o , 


o ft p o £ ? 

ffl o p 13 o o o 

C ffl ffl ffl p ft p 

P E w ft w p 

w 3 3 

•H o C c o ffl 

£ ffl o P ft 

W P ffl p -rt p 

ffl ft ffl W 13 

3 ffl Eh £ W ffl P 

O’ ft p w o , 

p P w E 3 

C • ffl w P 

ft e w 3 c ffl ffl 

o ffl p o' rt Pft- 

ffl ffl rt p p ffl o 

p £ o c p * E 

P hC ft 

C ffl OP - p 

o ft ft ffl w x rt ’ 

p o p rt P 3 

P ft P POP 

rt C rt p ft £ ft • 

P <D ffl rt ffl 

3 P w O P P O 

E ft ffl p p w C 

p P p rt rt O 

W Q) >> ft o 

ft O p ffl 

ft p p rt ft ft ft 

P «H C P P ffl 

o P P rt p N 

ft rt o <h ^ p 

ft ffl ft o p 

p p ft C ft rt * 

o w rt w <d P -i 


(ffl C O o be C 

: p E ft o p o 

> ft rt <d p ft rt 
p p p 

J rt bJD rt rt S w he I 

) 0 P P O P C 1 

) P P ffl 3 p p P i 

( rt ft P E <h 3 ft -t 


> ft P p p w > 

i ffl rt o rt rt ffl 0 ) 

ft ffl ffl p ft w x 
i 3 w p 3 <D 

IP (DP ffl P ft P 

! O p ft E o Eh O 


3 (DPP 
ffl w in Q) c 

I I <d ft £z> 

K W c E 

p in ffl ^ 

< f- w > p 

Ctf o o 

WWDSCS 


rt ft rt 
p • p ft 
, w P w 
U ffl rt ft 
ffl p ffl ffl 


p ft E p ‘ 

+j 0 W 

ft C rt 

COOP 


I O £ 5* 

ft 0 O 0 

I ffl p ft p 1 

W P, W «H 


^ W ffl P 

p w o . 
w E 3 



- ffl ft C 

w o 3 ffl l 

a o o wj c 

O P o 

pft rt 

rt M* M i 

P O p C I 

3 rrt P P I 
E P 3 ft -i 


rt rt ffl ffl 
P ft W x ■ 
C . ffl 
ffl p ft p 
E 0 Eh O 




